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Abstract 

Oscillations  of  clamped-clamped  and  free-free  micro-electromechanical  resonators  and 
resonator  arrays  have  been  studied  in  this  effort.  Piezoelectric  actuation  is  used  to  excite 
these  resonator  structures  on  the  input  side  and  piezoelectric  sensing  is  carried  out  on  the 
output  side.  Composite  structural  models  have  been  developed  for  these  filters,  and 
analyses  has  been  carried  out  to  explain  experimental  observations  of  nonlinear 
phenomena  as  well  as  to  guide  the  design  of  these  filters  are  presented.  Semi-analytical 
design  tools  for  micro-electromechanical  resonators  and  micro-electromechanical 
resonator  arrays  have  been  developed.  The  phenomenon  of  intrinsic  localized  modes  in 
resonator  arrays  has  also  been  studied,  and  it  is  shown  that  these  modes  can  be  explained 
as  forced  nonlinear  vibration  modes.  The  research  findings  can  open  the  doors  to  new 
resonator  array  designs. 


Keywords:  microscale  resonators,  buckling,  nonlinear  phenomena,  axial-load  effects, 
piezoelectric  actuation,  Duffing  oscillator,  resonator  arrays,  intrinsic  localized  modes. 

1.  Introduction 

Two  types  of  piezoelectrically  actuated  micro-scale  resonators,  which  are  attractive  for 
communication  and  signal-processing  applications  [1],  are  considered  in  this  research 
effort.  One  type  of  resonators  will  be  referred  to  as  the  AlGaAs  resonator,  and  the  other 
type  of  resonators  will  be  referred  to  as  the  PZT  resonator.  Both  types  of  resonators  are 
composite  structures,  and  the  PZT  resonators  have  asymmetric  cross-sections,  as  shown 
in  Figure  1.  During  the  fabrication  of  the  resonators,  residual  stresses  are  likely  to  be 


introduced  and  the  effect  of  these 
stresses  is  explored  in  this  work. 

As  pointed  out  in  the  authors’ 
recent  work,  the  resonators  also 
exhibit  non-linear  characteristics 
[2-5].  These  characteristics 
include  Duffing  oscillator  like 
response  during  resonance 
excitations  [6],  temporal  harmonics 
in  the  response,  and  spatial  patterns 
during  forced  oscillations  that 
cannot  be  explained  by 
conventional  linear  analysis. 

The  rest  of  this  report  is  organized 
as  follows.  In  Section  2,  the  first 
component  of  the  research  activity 
is  described.  This  section  contains 
a  discussion  of  a  semi-analytical 
tool  developed  to  study  transverse 
free  vibrations  of  clamped-clamped 
resonators  subjected  to  constant 
axial  loads.  In  Section  3,  nonlinear  analysis  used  to  study  dynamic  buckling  in  microscale 
resonators  is  presented.  Following  that,  progress  made  towards  using  an  experimentally 
obtained  frequency  response  for  parametrically  identifying  a  nonlinear  oscillator  model 
of  a  micro-resonator  is  described.  In  Section  5,  work  being  conducted  with  arrays  of 
composite  microresonators  is  briefly  discussed.  A  summary  is  provided  in  Section  6. 
References  including  the  publications  that  came  out  of  this  work  are  included  in  Section 
7.  The  research  participants  are  provided  in  Section  8  and  transitions  and  interactions  are 
discussed  in  Section  9.  A  list  of  publications  that  have  followed  from  this  work  is 
provided  in  Section  10,  and  invention  disclosure  made  through  this  work  is  provided  in 
Section  1 1 . 
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Figure  1 :  (a)  SEM  of  a  PZT  resonator  (Courtesy,  Maryland 
MEMS  Laboratory)  and  a  schematic  showing  the  details  [  1  ]. 


2.  Semi-Analytical  Tool  Development 


Based  on  a  formulation  with 
geometric  nonlinearities,  a  semi- 
analytical  tool  [7,  8]  has  been 
developed  for  studying  transverse 
free  vibrations  of  clamped-clamped 
and  free-free  resonators  subjected 
to  constant  axial  loads.  It  is  shown 
that  the  consideration  of  axial  loads 
is  important  to  predict  the  natural 
frequencies  of  the  resonators 
observed  in  the  experiments.  In 
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Figure  2:  Response  of  a  400  //m  AlGaAs  for  three  different 
values  of  the  axial  load  [7]. 


addition,  the  developed  model  and  the  numerical  implementation  can  be  used  to 
understand  the  influence  of  uncertainties  associated  with  the  fabrication  process. 
Representative  results  obtained  for  a  400  pm  long  resonator  are  shown  in  Figure  2. 
These  results  show  how  the  center  frequency  of  a  filter  can  be  shifted  by  using  an  axial 
load.  Additional  details  of  the  work  can  be  found  in  the  paper  [8]  included  at  the  end  of 
the  report. 

3.  Dynamic  Buckling  of  Composite  Micro-Scale  Structures:  Nonlinear 
Analysis 

The  resonators  considered  in  this  effort  are  based  on  the  piezoelectric  effect  (see  Figure 
1).  In  the  case  of  the  PZT  resonators,  the  elastic  substrate  is  a  SiC>2  layer,  on  top  of 
which  a  platinum  electrode  layer  is  deposited  throughout  the  length  of  the  structure.  A 
thin  layer  of  sol-gel  piezoelectric  film  is  located  on  the  top  of  this  electrode  layer.  To 
complete  the  structure,  another  platinum  layer  on  the  top  of  this  piezoelectric  film 
extends  over  one  quarter  of  the  length  from  each  anchor  and  the  mid-section  of  the 
resonator  structure  is  free  from  this  platinum  electrode  [1].  Due  to  the  asymmetry  of  the 
cross  section,  the  position  of  the  piezoelectric  layer  is  offset  from  the  neutral  axis,  and  in 
addition,  (tensional)  residual  stress  may  also  be  introduced  in  each  layer  during  the 
manufacturing  process.  In  some  typical  uses  of  this  resonator,  the  structure  is  driven  close 
to  its  first  resonance  frequency  with  a  DC  bias  in  the  input.  So,  the  axial  loads  in  a 
resonator  can  be  attributed  to  the  DC  bias  as  well  as  the  residual  stresses. 

It  is  believed  that  some  of  the  experimental  observations  made  in  previous  efforts  can  be 
explained  by  considering  oscillations  about  a  non-flat  equilibrium  position.  Here,  it  is 
assumed  that  this  non-flat  equilibrium  position  arises  due  to  buckling.  Hence,  the 
primary  thrust  of  this  component  has  been  to  examine  if  oscillations  of  a  buckled  system 
can  be  used  to  explain  experimentally  observed  nonlinear  motions.  Representative 
results  obtained  through  the  analysis  [9,  10,  11]  are  shown  in  Figure  3.  These  results  are 
in  good  agreement  with  the  corresponding  experimental  observations.  A  more  full 
discussion  of  the  analysis  is  provided  in  the  papers  [10,  11]  included  at  the  end  of  this 
report. 


mode,  (b)  Predicted  spatial  pattern  when  b  =  4. 65*  Iff4,  where  the  corresponding  natural  frequency  is  334 
kHz,  and  (c)  Analytical  prediction  for  the  frequency-response  curve  when  b  =  4.65  x  /  (T4,  p  =  4.65  x  10'1 , 
and  p  =  8.23x10“* ;  the  solid  line  represents  the  stable  branch  and  dashed  line  represents  the  unstable 
branch. 


4.  Experimental  System  Identification 


The  investigators  have  been  able  to  use 
experimentally  obtained  nonlinear 
frequency-response  data  to  determine  the 
linear  and  nonlinear  parameters 
governing  a  micro-resonator  [12,  13]. 
Representative  results  are  shown  in 
Figure  4.  Through  this  work,  the 
different  parameters  such  as  modal  mass, 
linear  stiffness,  nonlinear  stiffness, 
equivalent  damping  coefficient,  modal 
forcing,  and  residual  stress  are  identified 
for  a  microresonator  excited  close  to  its 
first  natural  frequency.  The  variations  of 
these  parameters  with  respect  to  the 
operating  conditions  have  also  been 
studied.  Full  details  of  the  identification 
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Figure  4:  Experimental  data  fit  compared  with  prediction 
from  forced  Duffing  model. 


scheme  are  included  in  the  paper  [13],  which  is  included  at  the  end  of  this  report. 


5.  Analyses  of  Microresonator  Arrays 

Both  linear  and  nonlinear  analyses  have  been  carried  out  to  determine  the  responses  of 
micro-resonator  arrays  to  different  excitation  conditions.  In  the  linear  analysis  [e.g.,  8], 
the  admittance  functions  of  clamped-clamped  AlGaAs  resonator  arrays  (e.g..  Figure  5)  is 
being  studied  to  aid  the  analysis  and  design  of  these  devices. 

In  the  nonlinear  analyses  [14],  the  authors  have  recently  examined  the  possibility  of 
realizing  intrinsic  localized  modes  (ILMs)  in  micro-resonator  arrays.  The  hypothesis  that 
these  modes  may  be  forced  nonlinear  vibration  modes  of  a  resonator  array  has  been 
analytically  and  numerically  studied,  and  it  has  been  shown  that  can  be  realized  as  forced 
nonlinear  vibration  modes.  Additional  details  on  this  analysis  is  included  in  the  paper 
[14]  included  with  this  report. 

6.  Summary 

It  is  believed  that  the  numerical  and  analytical  efforts  discussed  in  this  report  can  be  used 
as  a  basis  to  understand  nonlinear  phenomena  such  as  dynamic  buckling  and  nonlinear 
vibration  modes  in  composite  microscale  resonators  and  resonator  arrays,  parametrically 
identify  microresonators,  and  develop  semi-analytical  design  tools  for  microresonators 
and  microresonator  arrays.  The  present  work  can  help  take  advantage  of  nonlinearities 
for  designing  RF  filters,  mixers,  switches,  signal  processors,  and  other  MEMS  devices. 


Figure  5:  Illustrative  plots  of  the  admittance  functions  at  different  output  ports  of  a  microresonator  array 
with  respect  to  different  system  parameters. 
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Buckling  and  Free  Oscillations  of 
Composite  Microresonators 

He  Li  and  B.  Balachandran,  Fellow ;  ASME 


Abstract — Free  oscillations  of  piezoelectric,  microelectrome- 
chanical  resonators  are  considered  in  this  effort.  These  resonators 
are  modeled  as  clamped-clamped  composite  structures,  with 
stepwise  varying  properties  across  the  length  of  the  resonator.  The 
different  features  of  the  model  development  are  discussed,  and 
buckling  in  these  resonators  is  studied.  A  nonlinear  analysis  con* 
ducted  to  study  oscillations  about  a  buckled  position  is  presented. 
The  results  of  the  analysis  are  found  to  compare  well  with  the 
experimental  observations.  [1482] 

Index  Terms — Buckled  beam,  clamped-clamped  piezoelectric 
resonator  structure,  dynamic  buckling,  microelectromechanical 
systems  (MEMS)  resonator. 


I.  Introduction 

SCILLATIONS  of  microelectromechanical  resonators 
fabricated  as  clamped-clamped  composite  structures  are 
studied  here.  These  microresonators,  which  are  to  be  used 
as  filters,  are  important  for  mobile  communication  systems 
and  signal  processing  applications  (e.g.,  Fourier  transform 
computations)  [1].  Analytical  studies  of  the  nonlinear  dynamic 
response  of  electrostatically  actuated  microresonators  modeled 
as  beam-like  structures  have  been  recently  carried  out  [2].  The 
resonators  considered  in  this  effort  are  based  on  the  piezoelec¬ 
tric  effect,  as  shown  in  Fig.  1.  The  elastic  substrate  is  a  SiC>2 
layer,  and  on  the  top  of  this  layer,  a  platinum  electrode  layer  is 
deposited  throughout  the  length  of  the  structure.  A  thin  layer  of 
sol-gel  piezoelectric  film  is  located  on  the  top  of  this  electrode 
layer.  To  complete  the  structure,  the  top  layer  is  a  platinum 
electrode  layer  that  extends  over  one  quarter  of  the  length  from 
each  anchor  [3],  [4].  Each  resonator  structure  has  three  layers 
in  the  mid-span  where  the  top  electrode  layer  is  absent  and  four 
layers  elsewhere.  Due  to  the  asymmetry  of  the  cross  section, 
the  position  of  the  piezoelectric  layer  is  offset  from  the  neutral 
axis,  and  in  addition,  residual  stress  may  also  be  introduced  in 
each  layer  during  the  manufacturing  process. 

Here,  the  resonators  considered  typically  range  in  lengths 
from  1 00  /im  to  400  /xm,  with  a  width  of  20  /xm  and  a  total  thick¬ 
ness  of  about  2.3  /zm.  For  comparisons  of  model  predictions 
with  experimental  results,  particular  attention  is  paid  to  a  200 
/xm  long  resonator.  The  layer  thickness  values  for  this  resonator 
are  provided  in  Table  I.  As  shown  in  Fig.  2,  each  resonator  is 
modeled  as  a  composite  beam  with  stepwise  axially  varying 
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properties.  The  axial  stiffness  and  bending  stiffness  values  for 
a  resonator  considered  in  this  study  are  given  in  Table  II.  The 
subscripts  used  for  the  different  stiffness  values  correspond  to 
the  sections  shown  in  Fig.  2. 

II.  Experiments  and  Observations 

For  a  200-/zm-long  resonator,  the  first  natural  frequency,  was 
experimentally  determined  to  be  close  to  334  kHz,  while  the 
model  prediction  for  this  resonator,  without  consideration  of  the 
axial  force  and  nonflat  equilibrium  position,  is  about  1 86  kHz.  A 
sketch  of  the  experimental  arrangement  used  to  study  forced  os¬ 
cillations  of  this  resonator  is  shown  in  Fig.  3.  A  laser  vibrometer 
is  used  to  get  a  measure  of  the  transverse  vibrations  at  the  mid¬ 
point  of  the  resonator.  The  excitation  signal  fed  into  the  input 
port  of  the  resonator  consists  of  a  sinusoidal  component  close 
to  the  first  natural  frequency  and  a  dc  bias  offset.  In  Fig.  4,  a 
representative  spatial  response  distribution  of  this  resonator  ob¬ 
served  in  the  experiments  is  shown  [4].  This  spatial  distribution 
was  measured  by  using  a  scanning  laser  vibrometer.  The  exper¬ 
imental  observations  suggest  that  the  oscillations  may  be  taking 
place  around  a  nonflat  equilibrium  position. 

Here,  the  hypothesis  that  the  nonflat  equilibrium  position 
is  caused  by  buckling  is  explored  to  explain  the  experimental 
observation  shown  in  Fig.  4  and  accurately  predict  the  exper¬ 
imentally  obtained  value  of  the  first  natural  frequency.  This 
hypothesis  is  motivated  by  prior  work  conducted  with  buckled 
microstructures  [5]— [9]  and  large-scale  structures  [10],  [11]. 
These  prior  studies  on  electrostatically  actuated  microstructures 
have  by  and  large  focused  on  the  static  case,  and  in  addition, 
the  modeling  of  the  spatial  information  has  not  been  carried 
out  in  the  dynamic  case.  Furthermore,  in  all  of  the  previous 
studies,  the  structures  considered  have  a  uniform  section  unlike 
the  microstructures  considered  in  the  present  work.  Here,  the 
piezoelectrically  actuated  microstructures  have  stepwise  axially 
varying  properties  (see  Table  II). 

To  explain  the  response  characteristics  such  as  those  shown 
in  Fig.  4,  here,  the  buckling  of  a  beam  with  stepwise  axially 
varying  properties  is  studied  and  the  free  oscillations  about  dif¬ 
ferent  buckling  modes  are  examined.  The  model  development 
is  presented  in  Section  III,  along  with  the  buckling  analysis  and 
treatment  of  free  vibrations  about  a  buckled  position.  This  work 
is  of  general  nature  and  it  can  be  used  to  study  buckling  in  any 
composite  beam  with  stepwise  axially  varying  properties  and 
free  oscillations  of  such  a  structure  about  a  buckled  position. 
In  Section  IV,  the  model  predictions  for  microresonators  are 
provided,  compared  with  experimental  observations,  and  dis¬ 
cussed. 
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(a) 

sense 


(b) 

Fig.  1 .  Piezoelectric  microresonator:  (a)  SEM  picture  of  a  400-/im-long  PZT  resonator  (courtesy,  Maryland  MEMS  Laboratory)  and  (b)  schematic  of  piezoelectric 
resonator  [3]. 


TABLE  I 


Layer  Thickness  Values  for  a  200  /im  Composite  Resonator 

SiOj 

Bottom  Pt 

PZT 

Top  Pt 

(pm| 

|pm| 

|pm) 

||im) 

1.030 

0.085 

1.09 

0.090 

Fig.  2.  Clamped-c lamped  composite  beam  with  stepwise  axially  varying 
characteristics. 


TABLE  II 

Axial  Stiffness  and  Bending  Stiffness  Values  for  a 
200  /im  Long  Resonator 


EA, 

EA% 

EA> 

EI, 

Eh 

Eh 

|N| 

[N] 

IN) 

[Nm2] 

|Nm2] 

[Nm>| 

3.17 

2.88 

3.17 

1.39X10*12 

0.83X1012 

|  39X10  '2 

Spectrum  Analyzer  La mr  Doppler  VIbro meter 


Fig.  3.  Experimental  arrangement  showing  how  a  laser  vibrometer  is 
positioned  to  examine  transverse  vibrations  of  resonator.  The  resonator  is 
excited  by  signals  input  to  the  drive  electrode. 


Fig.  4.  Laser  vibrometer  measurement  of  a  spatial  pattern  observed  in 
experiments.  The  presence  of  spatial  harmonics  distorts  the  spatial  pattern  from 
the  typical  mode  shape  associated  with  the  fundamental  mode  of  vibration  of  a 
clamped-clamped  beam  [4]  (courtesy,  Maryland  MEMS  laboratory). 
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III.  Modeling  and  Quantitative  Analysis 
A.  Governing  Equations 

In  Fig.  2,  a  model  of  the  resonator  as  a  clamped-clamped 
beam  has  been  presented.  The  axially  varying  properties  along 
the  resonator  can  be  expressed  in  the  form 

n 

m(x)  =  ^™kHk{x), 

fc= l 
n 

EA{x)  =  YJEAkHk(x), 

k=  1 
n 

EI(x)=  J2EhHk(x)  (1) 

k=  1 

where  rrifc,  Elk,  and  EAk  are  constants,  x  is  the  spatial  variable, 
Hk(x)  =  [u(x  —  Xk-i)  -  u(x  —  x fc)]  is  the  Heaviside  function 
and  u(x)  is  the  unit  step  function,  and  n  is  the  total  number 
of  beam  sections,  which  is  3  here.  The  structure  extends  from 
x o  =  0  to  xn  =  /,  where  l  is  the  unbuckled  structure’s  length. 

For  the  considered  composite  structures,  taking  into  account 
that  the  length  to  thickness  ratio  is  more  than  ten,  the  width  is 
much  larger  than  the  thickness  in  each  section,  it  is  assumed 
that  each  section  of  the  structure  can  be  treated  as  an  Euler- 
Bemoulli  beam  and  the  lateral  displacement  along  the  y  di¬ 
rection  of  Fig.  1(b)  can  be  ignored.  Then,  for  initially  straight 
beams  undergoing  undamped  and  unforced  motions,  the  non¬ 
linear  equations  of  motion  can  be  written  as  [12] 


d2uk  (x,t)  d2uk(x,t) 

‘  EA"—~ 

k  =  1, 2, . . .  ,n  (2) 


d2wk(x,t)  -  d2wk{xyt)  d4wk(x,t) 

DP  ~  N°  W  +  Eh  di< 

=[“*-*•]  I ■ 

fc  =  1, 2, . . .  ,n  (3) 


where  the  nonlinear  axial  strain  has  the  form 


ek(x,t) 


k 


ds  —  dx 
ds 

duk  _  / du k\ 2  1  / du>fc\2 
dx  V  dx  )  2  V  dx  )  ’ 

1,2, ...,n. 


(4) 


In  (2)— (4),  the  symbol  Q  has  been  used  to  indicate  a  dimensional 
variable,  Wk(x,t)  is  the  transverse  displacement  of  the  beam 
in  the  kth  section,  Uk{x ,  t)  is  the  displacement  along  x  axis  in 
the  A;th  section,  t  is  the  time  variable,  and  N0  is  the  prescribed 
initial  static  axial  load,  which  is  assumed  to  be  a  constant  with 
respect  to  x.  In  (4),  dx  is  the  length  of  a  differential  element  of 
the  undeformed  beam  and  ds  is  the  corresponding  length  of  a 
differential  element  of  the  deformed  beam.  For  the  convenience 


of  further  analysis,  the  following  dimensionless  variables  are 
introduced: 


x  =  i' 


w 

w  =  v 


u=-r 


tk 


i  [EE 

l2  \  mk  ’ 


Nk  = 


N0l2 

EIk' 


(5) 


In  terms  of  the  dimensionless  spatial  variable  x,  the  boundaries 
of  the  structure  are  located  at  xo  =  0  and  xn  =  1.  By  using  the 
nondimensional  variables  given  in  (5),  (2)-(4),  can  be  written  in 
the  following  form: 

d2uk 


rl 


d2uk 

dt2 


rl 


d2wk  AT  d2wk 
_  dt2k  k  dx 2 


(6) 


=  [1  -  r2kNk\ 


d  \  dwk 
x^|efc- 


dx 


k  =  1,2, . . .  ,n. 


(7) 


e/k  (x,tk)  = 


duk 


_  (  duk  y 

dx  \  dx  ) 

1  ( dwk \ 2 

2  \dx  )  ’ 
k  =  1, 2, . . .  ,n. 


(8) 

In  (6)-(8),  the  parameter  rk  =  ( yj EIk/ EAk/l )  is  the  slender¬ 
ness  ratio  of  the  k\h  section  of  the  resonator. 

For  the  resonator  with  a  length  of  200  /im  and  the  stiffness 
values  provided  in  Table  II,  the  kth  section’s  slenderness  ratio 
rk  is  of  the  order  of  10~3;  this  means  that  the  longitudinal  inertia 
term  is  small  compared  to  the  longitudinal  stiffness  term  in  (6). 
Here,  based  on  the  small  value  of  the  slenderness  ratio,  it  is 
assumed  that  uk  =  0{w\)  [12].  Then,  (8)  can  be  reduced  to 

2 

+  ....  (9) 


_  duk  1  /  dwk  \ 

dx  2  \  dx  J 


ek(tk). 


(10) 


Dropping  the  terms  with  rk  and  the  cubic  terms,  (6)  can  be 
integrated  with  respect  to  the  spatial  variable  to  obtain 

duk  _  1  / dwk\2  _ 
dx  ~  2  V  dx  )  ' 

From  (10),  it  is  noted  that  the  displacement  uk  can  be 
regarded  as  a  quasi-static  motion,  and  that  the  expression  for 
ek  is  equivalent  to  the  expression  for  ek  in  (9),  which  has  the 
same  form  as  the  Green-Lagrange  strain  fn  [13].  Spatially 
integrating  (10)  one  more  time  and  evaluating  the  resulting 
expression  at  x  =  xkyek  can  be  determined  as 
uk  ( xk,tk )  -  uk  (xfc_i ytk) 


£k(tk)  = 


2  (xfc 

k  =  1,2, ...  ,n. 


(Xk  -  Xfc-l) 

j _ r 

-Xk-l)  Jx k 


(11) 
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From  ( 1 0),  it  had  been  noted  that  ek  is  constant  with  respect  to 
x.  Based  on  this  observation,  the  force  equilibrium  condition  at 
the  interconnection  of  any  two  adjacent  sections  requires  that  the 
axial  load  due  to  this  strain  is  constant  among  all  of  the  sections; 
that  is 


Q{t)  =  EAk£k{i ).  (12) 


After  substituting  (11)  for  ek  into  (7)  and  simplifying,  the 
resulting  equation  of  motion  governing  the  transverse  vibrations 
of  the  microresonator  with  stepwise  axially  varying  properties 
is  obtained  as 


d2Wk  (  d4Wk 

~dq  ' 


dx 4 


Nk(tk)  + 


Uk(Xkytk)  -  Uk(Xk-l,tk) 


2r2k(xk 

k  =  1,2, 


rUxk  ~  Zfc-l) 

I _  f 

-  Xk-l)  Jx 


(£)' 


dx 


CpWk 

dx2 


=  0, 
(13) 


For  the  particular  case  of  a  clamped-clamped  resonator,  the 
associated  boundary  conditions  are 

w(x,  t)  =  0  and  )  —  0  at  x  =  0  and  x  =  1.  (14) 

ox 

The  governing  equation  given  by  (13)  is  different  from 
that  treated  in  previous  work  [10],  [11],  since  the  resonator 
has  axially  varying  properties.  This  equation  along  with  the 
boundary  conditions  represents  a  nonlinear  model  that  can  be 
used  to  study  transverse  vibrations  of  a  beam  with  stepwise 
axially  varying  properties  and  sections  that  have  “small”  slen¬ 
derness-ratio  values. 


B.  Static  Buckling  Problem 


To  solve  (13),  the  static  buckling  problem  is  first  addressed. 
It  is  assumed  that  the  dimensionless  static  buckling  deflection 
of  the  resonator  can  be  expressed  as 

n 

X(x)  =  ^2xk(x)Hk{x) 

k=l 

Tl 

7T(x)  =  k(x)Hk(x)  (15) 

k=l 


where  Xfc(x)  and  ^k(x)  are  the  transverse  and  axial  displace¬ 
ment  of  the  A:th  section  of  the  beam,  respectively.  For  the  static 
buckling  analysis,  the  inertia  term  in  (13)  is  dropped,  and  the 
static  compressive  axial  forces  —  Pq  and  —Po,k  are  substituted 
for  No  and  Nk  respectively  in  (5)  and  (13).  After  carrying  these 
substitutions  through,  the  nonlinear  equation  governing  the 
equilibrium  of  the  beam  subjected  to  a  static  axial  force  can  be 
written  as 


d4Xk 
dx 4 


+ 


p  _  TTfcfofc)  ~  frfcQgfc- l) 

rk(Xk  ~  xk- 1) 

1 


2 r£(arfc  -  xk- 


dx 


d?Xk 
dx 2 


=  0,  k  =  1, 2, . . .  ,n. 


Hx,i)  :  Dynamic  displacement  Total  dltplacwnant 


The  associated  boundary  conditions  follow  from  (14)  as 

X(x)  =  0  and  ^  ^  =  0  at  x  =  0  and  x  =  1.  (17) 

dx 

In  order  to  obtain  the  solution  of  (16),  the  critical  buckling 
force  must  be  first  calculated  by  attacking  the  linear  buckling 
problem  for  this  composite  resonator  with  axially  varying  prop¬ 
erties.  To  this  end,  after  switching  to  the  dimensional  variable 
x ,  the  critical  buckling  configuration  ^(x)  is  determined  with 
respect  to  the  applied  critical  axial  load  Pc.  First,  it  is  assumed 
that 

n 

V'(i)  =  '%2'Pk(x)Hk(x).  (18) 

k=l 


Then,  noting  that  each  of  the  k  sections  has  uniform  proper¬ 
ties,  the  linear  static  buckling  problem  is  written  as 

+  =  °’  fc  =  1’2,..,n-  (19) 

Introducing  the  dimensionless  variables  from  (5)  into  (19) 
results  in 

d2V»fc  _n  ,  _  ,  O 

dx4  +<kdx2  ~0,  fc-1-2,...,n  (20) 

where  =  PCyk  —  (PJ2 / Eh)-  The  general  solution  of  (20) 
is  given  by 

V'fc(x)  =  aiik  +  a2ykx  +  a3ik  cos(Ckx)  +  a4fk  sin((kx).  (21) 


To  determine  the  arbitrary  constants  ajyk  for  j  =  1,2,3, 4 
and  k  =  1, 2, . . . ,  n,  the  following  boundary  conditions  for  a 
clamped-clamped  beam  are  used 

.  v  dib\(x} 

\j)\ ( x )  =  0  and  — - - =  0  at  x  =  xq, 

dx 

djl )  (x) 

\pn(x)  =  0  and  — j - =  0  at  x  =  xn.  (22) 

dx 

In  addition,  as  the  beam  is  a  continuous  system,  two  adja¬ 
cent  sections  have  to  share  the  same  transverse  displacement, 
same  slope,  same  bending  moment,  and  same  shear  force  at  the 
respective  connection  point.  This  leads  to  the  following  nondi- 
mensional  form  of  the  compatibility  conditions  [14] 


V'fc-i(ar)  =tpk(x), 

dipk- i(ar)  _  drpkjx) 
dx  dx 


EIk- 1 


dVfc-i(x) 

dx2 


-  Elk 


d2ipk{x) 
dx2  ’ 


EIk  <Pj>k-l(x)  _EIkf^k{x) 


dx 3 


dx 3 


at  x  =  xk_\ ,  k  =  2, 3, . . .  ,n. 


(23) 


x 


06) 
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Fig.  6.  (a)  First  static  buckling  mode  shape  of  the  200  /im  PZT  resonator,  (b)  first  natural  frequency  versus  6,  (c)  first  mode  shape  of  free  vibration  about  first 

buckling  mode  for  6  =  0.001,  and  (d)  first  mode  shape  of  free  vibration  about  first  buckling  mode  for  b  =  0.02. 


For  n  =  3,  after  substituting  (21)  into  (22)  and  (23),  one  ob¬ 
tains  a  set  of  twelve  homogeneous  algebraic  equations  that  de¬ 
fines  an  eigenvalue  problem  for  the  ajyk  and  the  eigenvalues  (fe¬ 
lt  is  important  to  note  here  from  the  definition  of  (fe  provided  in 
the  context  of  (20),  that  the  different  eigenvalues  are  functions 
of  the  critical  buckling  force.  From  (9)  and  (12),  it  is  noted  that 
in  the  kih  section,  the  axial  displacement  #fc,  the  transverse  dis¬ 
placement  rpi t  and  the  longitudinal  extension  measure  eCyk  have 
to  satisfy 


post-buckling  deformation  with  respect  to  the  shape  of  the  crit¬ 
ical  buckling  mode.  From  (9)  and  (24)  along  with  (26),  one  can 
obtain  that 

wfc  -  b\tik 

ek  =  b\eCtk.  (27) 

Next  substitution  of  (20),  (24),  and  (26)  into  (16)  results  in 


dtik  1  /  dipk  \  2 

(24) 

Cc  fc "  dx  +  2  (  dx  ) 

Qc(t)  =EAkeCtk(t). 

(25) 

After  calculating  the  critical  buckling  force,  the  post-buck¬ 
ling  problem  is  considered.  When  the  compressive  axial  force 
is  larger  than  the  critical  buckling  force,  the  linear  buckling 
problem  given  by  (20)  cannot  be  used  to  study  the  beam’s  defor¬ 
mation.  The  nonlinear  equation  given  by  (16)  needs  to  be  con¬ 
sidered. 

For  solving  (16),  it  is  assumed  that 


6 


2 

k 


rljPoj,  -  PClk) 

&C,k 


k  =  1,2,.  ...n. 


(28) 


After  substituting  (25)  into  (28),  and  returning  to  dimensional 
variables,  it  is  found  that  6fe  is  the  same  constant  with  respect  to 
x  for  all  the  sections;  that  is 


b2  _  (Po  ~  Pc) 
Qc 


(29) 


Therefore,  the  solution  of  (16)  and  (27)  becomes 


Xk(x)  =  bkipk(x),  k  =  1, 2, . . .  ,n  (26) 

where  bk  is  a  nondimensional  factor,  which  is  called  the  buck¬ 
ling  factor.  This  factor  is  a  measure  of  the  amplitude  of  the 


Xk(x)  =brpk(x), 
nk(x)  -b2dk(x), 

ek=b2eCtk,  k  =  1, 2, . . .  ,n. 


(30) 
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Fig.  7.  (a)  The  third  static  buckling  mode  shape  of  the  200  /im  PZT  resonator,  (b)  first  natural  frequency  versus  6,  (c)  first  mode  shape  of  free  oscillation  about 

third  buckling  mode  for  b  =  0.001,  (d)  first  mode  shape  of  free  oscillation  about  third  buckling  mode  for  6  =  0.004,  and  (e)  first  mode  shape  of  free  oscillation 
about  third  buckling  mode  for  b  =  0.01. 


C.  Free  Linear  Vibrations  About  the  Postbuckled  Position 

Free  oscillations  of  the  undamped  resonator  modeled  as  a 
beam  with  stepwise  axially  varying  properties  are  considered 
in  this  section.  The  equation  of  motion  and  the  corresponding 
boundary  conditions  for  an  undamped,  unforced  beam  subjected 
to  an  axial  load  are  given  by  (13)  and  (14).  The  solution  of  this 
system  can  be  written  as  the  sum  of  a  static  component  and  a 
dynamic  component  [5]  (see  Fig.  5) 

Wk(x,  tk)  =  tyk(x)  +  Vk(x}  tk ), 

uk{x,tk)  =  b2dk(x)  +  AOMfc),  k  =  1,2 (31) 


where  Vk(x,  tk)  and  fk{x ,  tk)  are  the  A;th  component  of  the  dy¬ 
namic  deflection,  ipk{x)  and  flk{x)  are  the  kih  component  of  the 
static  buckling  mode  shape  ip(x)  and  il)(x),  respectively,  and  b 
is  the  buckling  level  factor  defined  by  (29). 

Next,  the  natural  frequencies  Cj  and  mode  shapes  associated 
with  free  oscillations  around  the  postbuckled  position  b^(x) 
are  examined.  To  find  the  natural  frequencies  and  mode  shapes 
from  (13),  (31)  is  substituted  into  (13).  It  is  assumed  that 
\vk{x,tk)\  <C  \bipk(x)\y  fk  =  0(v%),  and  the  separation  of 
variables  is  used  to  express  the  solution  as 

Vk(x,tk)  -  <pk{x)e'“ktl' 


(32) 
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Fig.  8.  (a)  The  fifth  static  buckling  mode  shape  of  the  200  /an  PZT  resonator,  (b)  first  natural  frequency  versus  b,  (c)  first  mode  shape  of  free  oscillation  about 

fifth  buckling  mode  for  b  =  0.001,  and  (d)  first  mode  shape  of  free  oscillation  about  fifth  buckling  mode  for  b  =  0.005. 


where  <t>k(x)  is  the  mode  shape  and  u>k  =  y/(mk/ EIk)l2v 
is  the  nondimensional  quantity  associated  with  section  k  = 
1,2, ...,n  Making  use  of  (20),  (13),  and  (32),  the  following 
equation  is  obtained: 


—  U)\4>k  +  P z,k 


d24>k  ,  d4<j>k 


dx 2 
EAkl2b2 


dxA 

d2xi)k 


EIk(xk  -xk-\)  dx2 
k  =  1,2, ...  ,n. 


rXk  r 

[dx  dx  \ 


(33) 


where 


(PcM  +  yJi*k+  4^2), 

A,k  =  \j 2  (--Pc,*:  +  \J Plk  +  4w£)  .  (36) 


After  substituting  (34)  and  (35)  into  (33),  the  differential 
equation  governing  the  particular  solution  (t>p,k(x )  is  obtained 
as 


The  general  solution  of  (33)  is  composed  of  a  homogeneous 
solution  <))tx,k{x)  and  a  particular  solution  4>Ptk(x);  that  is. 


4>k(x)  =  d>h,k(x)  +  <t>PAx)-  (34) 


As  in  the  critical  buckling  case,  the  spatial  displacement  term 
0fc(x)  has  to  satisfy  the  boundary  conditions  and  the  compati¬ 
bility  conditions  that  are  similar  to  those  given  by  (22)  and  (23). 
The  homogeneous  solution  is  given  by 


2i  .  d4<t>P,k  .  D  d2<t>p,k 

-  uk4>p,k  +  —jzr~  + 


dx 4 
EAkl2b2 


dx2 
d2ipk  rXi 


EIk(xk  -  xk- 1)  dx2 

EAkP fc2 

EIk(xk  -  Xk-i)  dx2 


J 

‘'H-l 


difrk  d(f>Pyk 


r 

JXk-i 


dx  dx 
d'lpk  d(f>h,k 


dx 


dx  dx 


dx. 


(37) 


It  follows  that  the  solution  of  (37)  can  be  written  as 


<t>h,k{x)  =  Ci.tsinfAi.fcx)  +  C2,fcCos(Ai,fcx) 

+C3,fc  sinh(A2,fcx)  +  C4,fc  cosh(A2-fcx)  (35)  <f>P,k(x)  =  C5,fc[a3,fc  cos((fcx)  +  a4,k  sin(Ctx)]  (38) 
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Fig.  9.  Predicted  spatial  pattern  for  b  =  0.00235,  where  the  corresponding  natural  frequency  is  333  kHz. 


where  the  k  are  constants.  Making  use  of  (37)  and  (38),  the 
relations  between  the  constants  Cy,*,  j  =  1, 2, . . . ,  5  are  ob¬ 
tained  as 

Si.fcCi.fc  +  B2,kC2,k  +  B$,kC$yk 

+#4,fcC4,fc  4-  BstkCstk  =  0  (39) 

where  the  Bj s  are  defined  in  Appendix  I.  Equation  (39)  con¬ 
stitutes  an  eigenvalue  problem  for  the  Cj and  the  natural  fre¬ 
quencies  l Ok-  After  determining  these  values  for  Cjtk>  the  ex¬ 
pression  for  the  free  vibration  mode  shape  of  the  postbuckled 
beam  with  stepwise  axially  varying  properties  can  be  written  as 


<f>{x)  =  y>t(x)fffc(j).  (40) 

fc=l 


IV.  Results  and  Discussion 

For  the  200-/xm-long  resonator  discussed  in  Sections  I  and 
II,  the  different  buckling  modes  determined  for  different  levels 
of  the  b  factor  are  shown  in  Figs.  6-9.  The  shapes  are  similar  to 
those  determined  at  the  macroscale  for  other  structures  [10]. 

As  discussed  in  Section  II,  the  first  natural  frequency  of  a 
200-/xm-long  resonator  was  determined  as  334  kHz  [4].  The 
amplitude  of  the  forced  oscillation  was  determined  to  be  in  the 
range  of  100  nm  to  1  fim  depending  on  the  drive  voltage  am¬ 
plitude  of  the  sinusoidal  signal.  In  Fig.  8(b),  the  variation  of  the 
first  natural  frequency  versus  the  scaling  parameter  b  is  shown 


for  free  oscillations  about  the  fifth  buckling  mode  of  the  res¬ 
onator.  When  b  is  equal  to  0.002  35,  the  first  natural  frequency 
of  the  resonator  is  determined  from  the  model  as  333  kHz,  which 
is  close  to  the  experimental  result;  it  is  important  to  note  that  for 
the  same  frequency  value,  the  corresponding  value  of  b  depends 
on  the  beam  stiffness  and  geometric  properties  as  well  as  how 
one  scales  the  static  critical  buckling  mode  shape.  In  this  work, 
the  buckling  factor  6  is  obtained  by  normalizing  the  nondimen- 
sional  static  critical  buckling  mode  shape  \p,  so  that  its  inner 
product  over  the  nondimensional  variable  running  from  0  to  1 
is  1.  The  corresponding  deflection  w  =  (hip  +  <p)  is  shown  in 
Fig.  9.  In  this  figure,  the  amplitude  value  of  (p  is  tuned  to  be 
10“3,  so  that  the  highest  displacement  over  the  spatial  coordi¬ 
nate  is  200  nm.  This  free  oscillation  mode  shape  is  similar  to  the 
spatial  patterns  observed  during  forced  oscillations  when  the  ex¬ 
citation  frequency  is  close  to  the  first  natural  frequency  of  the 
resonator  (see  Fig.  4). 

The  agreement  between  the  analytical  prediction  and  exper¬ 
imental  data  suggests  that  the  hypothesis  that  the  nonflat  equi¬ 
librium  position  of  the  resonator  is  caused  by  buckling  can  be 
a  valid  one.  Along  with  the  work  reported  in  reference  [5],  the 
present  study  provides  the  first  evidence  of  such  a  phenomenon 
in  microscale  resonators.  In  addition,  the  present  work  can  be 
used  as  a  basis  to  study  buckling  and  free  oscillations  of  res¬ 
onators  with  stepwise  axially  varying  properties. 


Appendix 

Here,  the  details  of  the  constants  that  arise  in  the  context  of 
(40)  are  provided.  See  the  equation  at  the  top  of  the  next  page. 
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B  i,fc  = 


Eh(xk  -  Xk-i) 

x  <  a2jfesin(A1,fca:)  +  a3,fcC*:Ai,fc 


-*  K  K 

1 02, k  sir 


cos[(Cit  -  Ai,fc)x]  cos[«fc  +  Ai,fc)x] 

2(0 k-Ai,fc)  2(Cfc  +  A1>fc) 


B2,k  = 


+<J4,fcCfc^l,fc 

EAkl2b2Pc,k 


sin[(Q  -  Ai,fc)x]  sin[(Cfc  +  Ai,fc)x] 
2  ((fc  —  Ai,it)  20  4-  Ai/t 


EIk(xk  -  xk-i) 

x  <  a2.fc  cos(Aitfcx)  +  a3tkCkh.k 


|a2,fc 


sin[(Q  -  Ai,fc)x]  _  sin[(Q  +  Ai,fc)x] 

2(0  -  Ai,fc)  2(0  +  A ,,*) 


Bz,k  = 


+a4,fcCfcAl,fc 

EAkWl^k 

EIk(xk  —  Xk-i) 


-  cos[(0  -  Ai,fc)x]  cos[(0  +  Ai,fc)x] 


2(0  -  Ai,*) 


2(0  +  Ai,/t) 


„  I  .  w,  V  ®3  j.(fcA2  jfc 

x  <  a2,fc  sinh(A2ifcx)  -  2  2 

{  a2  ,k  +  O 


A2,fc  sin(O^)  sinh(A2,fcx)-0  cos(Cfcx)  cosh(A2,fcx) 


+ 


#4,fc  — 


fl4,fcOA2,J; 

A^.fc  +  Cfc  . 

EAkWj^k 

EIk(xk  -  Xfc-i) 


A2,fc  cos(Ox)  sinh(A2,fcx)  +0  sin(Ox)  cosh(A2,*x) 


x  a3.fcOA2.fc 

1  Q’2,k  COSn(A2(fcl)  .  2  ,  y-2 

[  A2  ,fc  +  Sfc 


A2,fc  sin(O^)  cosh(A2,*;x)-0  cos(Ox)  sinh(A2>fcx) 


+ 


g4,fcQA2,fc 

A2,fc  +  Cfc  . 


A2,it  cos(Cfcx)  cosh(A  2,*x)  +0  sin(Cfcx)  sinh(A2jfcx) 


Bs.fc  = 


2^  EAkl2b2Pc  k 

~UJk  + 


EIk{xk  -  xk- 1) 

+  a2,*a3>fc  cos(C*x)  +  a2,fca4,fc  sin(Otx) 


3,<:  .  4’*Cfcsin(2Cfcx)  +  a3,fca4,fc  0  cos(20 x) 
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Abstract 

Free  and  forced  oscillations  of  piezoelectric,  microelectromechanical 
resonators  fabricated  as  clamped-clamped  composite  structures  are  studied 
in  this  effort.  Piezoelectric  actuation  is  used  to  excite  these  structures  on  the 
input  side  and  piezoelectric  sensing  is  carried  out  on  the  output  side.  A 
refined  integro-partial  differential  model  is  developed  for  a 
clamped-clamped  composite  beam  structure  and  used  for  studying  the 
nonlinear  transverse  vibrations  of  these  resonators.  This  model  accounts  for 
the  longitudinal  extension  due  to  transverse  vibrations,  distributed  actuation 
and  axially  varying  properties  across  the  length  of  the  structure.  Free 
oscillations  about  a  post-buckled  position  are  studied,  and  for  weak  damping 
and  weak  forcing,  the  method  of  multiple  scales  is  used  to  obtain  an 
approximate  solution  for  the  response  to  a  harmonic  forcing.  Analytical 
predictions  are  also  compared  with  experimental  observations.  The  model 
development  and  the  analysis  can  serve  as  a  basis  for  analysing  the 
responses  of  other  composite  microresonators. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Nonlinear  oscillations  of  microelectromechanical  resonators 
fabricated  as  clamped-clamped  composite  structures  are 
studied  in  this  paper.  These  microresonators,  which  are  to 
be  used  as  filters,  are  important  for  mobile  communication 
systems  and  signal  processing  applications  (e.g.,  Fourier 
transform  computations)  [1].  While  the  nonlinear  response 
of  a  microresonator  subjected  to  an  electrostatic  actuation  has 
received  considerable  attention  [2],  the  nonlinear  response  of 
a  microresonator  subjected  to  a  piezoelectric  actuation  has 
received  limited  attention  recently  (3]. 

The  resonators  considered  in  this  effort  are  based  on  the 
piezoelectric  effect,  as  shown  in  figure  1.  Fabrication  details 
for  these  resonators  can  be  found  in  [4].  The  dimensions  of  the 
resonators  considered  in  this  study  typically  range  in  length 
from  50  /zm  to  400  /zm,  with  a  width  of  20  /zm,  and  a  thickness 
of  about  2.3  /zm.  The  elastic  substrate  is  an  amorphous  SiC>2 
layer,  on  the  top  of  which  a  thin  platinum  electrode  layer  is 
deposited  first,  followed  by  a  layer  of  sol-gel  piezoelectric  film 
throughout  the  structure’s  length.  To  complete  the  structure,  a 
platinum  electrode  layer  that  extends  over  one  quarter  of  the 
length  from  each  boundary  is  deposited  as  the  top  layer  [5], 


so  that  each  resonator  structure  has  three  layers  in  the  mid¬ 
span  where  the  top  electrode  layer  is  absent  and  four  layers 
elsewhere.  Due  to  the  asymmetry  of  the  cross  section,  the 
position  of  the  piezoelectric  layer  is  offset  from  the  neutral 
axis,  and  in  addition,  residual  stress  may  also  be  introduced  in 
each  layer  during  the  fabrication  process. 

For  comparisons  of  model  predictions  with  experimental 
results,  particular  attention  is  paid  to  a  200  /zm  long  resonator. 
The  material  properties  and  thickness  values  for  the  different 
layers  are  provided  in  tables  1  and  2.  As  discussed  later  in  the 
third  section,  each  resonator  is  modeled  as  a  composite  beam 
with  properties  that  vary  in  a  stepwise  manner  from  section  to 
section,  as  shown  in  figure  2.  The  values  of  axial  stiffness, 
bending  stiffness  and  mass  per  unit  length  for  the  different 
sections  are  given  in  table  3  and  these  values  are  identified  by 
using  the  subscripts  used  in  figure  2. 

The  experimental  observations  that  will  be  used  for 
comparison  with  the  analytical  predictions  are  collected 
together  in  the  next  section.  In  the  third  section,  the 
model  development  carried  out  to  study  buckling  in  a 
composite  structure  is  presented,  and  in  the  fourth  section, 
analysis  of  linear  free  oscillations  of  the  structure  is 
detailed.  Forced  oscillations  are  examined  in  the  fifth  section 
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Figure  1.  Piezoelectric  microresonator:  (a)  SEM  of  a  50  /zm  long 
PZT  microresonator  [6]  and  (b)  a  schematic  showing  the  details  [5]. 


Table  1.  Material  properties  for  different  layers  of  a  200  /tm 
composite  resonator. 


Si02 

Bottom  Pt 

PZT 

Top  Pt 

P  (kg  m  3) 
E(  Pa) 

2194 

100  x  10” 

17  839 

160  x  I09 

8800 

25  x  109 

18  762 

160  x  109 

Table  2.  Thickness  values  for  a  200  /zm  composite  resonator. 


Si02  (/zm) 

Bottom  Pt  (/zm) 

PZT  (/zm) 

Top  Pt  (/zm) 

1.030 

0.085 

1.09 

0.090 

and  the  results  and  discussion  are  presented  in  the  sixth 
section. 


Figure  3.  Experimental  arrangement  of  laser  vibrometer  for 
examining  transverse  vibrations  at  the  mid-span  of  the  resonator. 
This  structure  is  excited  by  signals  input  to  the  drive  electrode. 


Table  3.  Values  of  axial  stiffness,  bending  stiffness  and  mass  per 
unit  length  for  the  200  /z m  long  resonator  of  interest. 


EAX  (N) 

EA2  (N) 

EA>  (N) 

3.17 

2.88 

3.17 

El,  (N  m2) 

Eh  (N  m2) 

Eh  (N  m2) 

1.39  x  10-'2 

0.83  x  10-'2 

1.39  x  10-'2 

m,  (kgm"‘) 

m2  (kg  m_l) 

m3  (kg  m_l) 

3.01  x  10-' 

2.68  x  10-' 

3.01  x  10-' 

2.  Experiments  and  observations 

In  figure  3,  a  sketch  of  the  experimental  arrangement  used 
to  study  transverse  vibrations  of  a  microresonator  is  shown. 
The  silicon  wafer  containing  the  resonator  is  first  placed  on  a 
probe  station,  and  the  drive  port  of  the  resonator  is  electrically 
connected  to  provide  a  sinusoidal  excitation  along  with  a  dc 
bias  input.  A  laser  vibrometer  is  used  to  measure  the  transverse 
vibrations  at  the  mid  point  of  the  resonator,  and  the  data  are 
collected  by  using  a  dynamic  spectrum  analyser. 

For  a  200  /zm  long  resonator,  the  first  natural  frequency 
was  experimentally  determined  to  be  close  to  313  kHz,  while 
the  model  prediction  for  this  resonator,  without  consideration 
of  the  axial  force  and  possible  nonflatness  of  static  state,  is 
about  1 86  kHz.  In  figure  4(a),  a  representative  spatial  response 
distribution  of  another  resonator  obtained  by  using  a  scanning 
laser  vibrometer  is  shown.  The  resonator  is  excited  by  a 


Figure  2.  Cl  am  ped-c  lamped  composite  structure  model  of  piezoelectrically  actuated  resonator  shown  in  figure  1 
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Figure  4.  Microresonator  response:  (a)  laser  vibrometer  measurement  of  a  spatial  pattern  observed  in  experiments;  ( b )  hardening 
Diiffing-type  frequency  response  of  a  200  /un  PZT  resonator  excited  by  a  sinusoidal  input  signal  with  an  amplitude  of  0.398  V  and 
(c)  spectrum  of  the  laser  vibrometer  measurement  at  the  mid-span  of  a  200  resonator,  in  response  to  an  excitation  with  a  frequency  close 
to  the  first  natural  frequency  and  an  amplitude  of  60  mV. 


sinusoidal  excitation  with  a  frequency  close  to  the  first  natural 
frequency  [6].  The  presence  of  spatial  harmonics  distorts  the 
spatial  pattern  from  the  typical  mode  shape  associated  with  the 
fundamental  mode  of  vibration  of  a  clamped-clamped  beam. 
A  typical  frequency  response  observed  in  the  experiments  is 
shown  in  figure  4 {b).  The  response  amplitude  varies  from 
a  few  hundred  nanometers  to  a  micrometer  depending  on  the 
excitation  level.  In  figure  4(c),  it  is  illustrated  that  the  response 
spectrum  consists  of  second  and  third  harmonics. 

The  presence  of  second  harmonics  suggests  that  the 
oscillations  may  be  taking  place  about  a  nonflat  static 
equilibrium  position  (e.g.,  [7]). 


In  this  paper,  the  hypothesis  that  the  nonflat  static 
equilibrium  position  is  caused  by  buckling  is  pursued  along 
the  lines  of  our  recent  work  ([3]).  This  hypothesis  is  motivated 
by  prior  work  conducted  with  buckled  microstructures 
[8-14]  and  large-scale  structures  [15].  The  prior  studies 
on  microstructures  have  by  and  large  focused  on  the  static 
cases  of  electrostatically  actuated  microstructures  [11-14], 
and  in  addition,  these  structures  were  modeled  as  uniform 
beams.  By  contrast,  in  this  work,  the  microstructures 
are  piezoelectrically  actuated  and  these  structures  have 
properties  that  vary  stepwise  in  the  axial  direction.  Here,  a 
refined  model  of  a  composite  structure  with  axially  varying 
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Figure  5.  Free-body  diagram  of  an  infinitesimal  beam  section.  The  beam  is  deformed  from  an  initial  flat  position.  Ax  is  the  length  of  the 
undeformed  beam  segment  and  As  is  the  length  of  the  corresponding  deformed  segment. 


properties  is  developed  and  used  to  predict  the  spatial 
responses  and  temporal  responses  of  microresonators  and  also 
accurately  predict  the  first  natural  frequency  of  the  considered 
microresonator. 

3.  Model  development  and  buckling  analysis 

Following  that,  this  model  is  extended  to  apply  to  a  structure 
with  stepwise  variation  of  the  properties  along  the  length  of 
the  structure,  such  as  that  shown  in  figure  2. 

This  model  is  used  to  examine  static  buckling  in  composite 
structures  such  as  the  microresonators  considered  in  this 
paper,  and  after  establishing  the  static-equilibrium  position 
as  a  buckled  position,  free  oscillations  are  considered  about 
a  buckled  position  in  the  fourth  section.  This  work  is  of  a 
general  nature  and  it  can  be  used  to  study  buckling  in  any 
composite  beam  with  axial  properties  that  vary  in  a  stepwise 
fashion  across  the  structure’s  length.  The  treatment  presented 
in  this  paper  is  a  refinement  of  an  earlier  treatment  presented 
by  the  authors  [17].  Compared  to  the  authors’  earlier  work, 
the  treatment  of  kinematics  is  different  in  this  case. 

3.1.  Governing  equations 

To  establish  the  governing  equations,  apart  from  the  Euler- 
Bemoulli  beam  assumptions,  it  is  also  assumed  that  the  axial 
and  transverse  displacements  are  small  compared  to  the  length 
and  that  the  initial  flat  position  shown  in  figure  5  is  free  of 
stresses  and  external  forces.  The  governing  equations  along 
axial  and  transverse  directions  can  be  obtained  as  [  19] 


where  the  carat  symbol  4A*  has  been  used  to  indicate  a 
dimensional  variable;  U(x,t)  and  W(x,t)  are  the  axial  and 
transverse  displacements;  x  is  the  spatial  variable  and  t  is  the 
time;  pA  is  the  mass  density  per  unit  length;  EA  is  the  axial 
stiffness;  El  is  the  bending  stiffness;  [lu  and  p.w  are  the  viscous 
damping  factors  for  motions  along  the  axial  and  transverse 
directions,  respectively.  As  shown  in  figure  5,  f,(x,t )  and 
/„( x,  t )  are  distributed  forces  along  the  tangent  and  normal 
directions,  and  m(Jt,  t)  is  a  distributed  moment;  the  axial  force 
N(x,  t )  is  given  by 


N (Jr,  t)  =  EAe,  (3) 

where  e  is  the  Green-Lagrange  strain  measure  [20], 
corresponding  to  the  elongation  of  the  beam  along  the 
longitudinal  direction.  The  expression  for  e  will  be  given 
later  in  this  paper. 
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a/2  a<  r2  s*2  a*4  a* 

Here,  noting  that  the  resonator  is  actuated  by  a  distributed 
piezoelectric  layer,  as  shown  in  figures  2  and  5,  the  boundary 
conditions  for  axial  and  transverse  vibrations  are  considered 
as  shown  in  figure  6. 

dW(x,t) 


W(x,  t)  =  0, 


W(xtt)  =  0, 


dx 

dW(xj) 

fix 


=  0,  U(x,  t)  =  0  at  jc  =  0, 


=  0,  U (jc,  t)  =  e  at  x  =  1. 


(7) 

Equations  (5)  and  (6)  along  with  the  boundary  conditions 
(7)  represent  the  governing  equations  of  a  clamped-clamped 
structure.  From  the  form  of  equation  (5),  it  is  clear  that 
the  longitudinal  dynamics  is  neglected  but  the  dependence 
of  the  longitudinal  strain  on  the  transverse  excitation  and 
displacement  are  captured. 

3.2.  Governing  equations  for  microresonators  with 
stepwise  varying  properties 

For  the  microresonators  shown  in  figure  2,  as  discussed  in  the 
previous  work  of  the  authors  [17],  the  axial  properties  can  be 


pA(x)  =  £/>/»*//*(*),  EA(x)  = 

*=l  k=\ 


(8) 


Figure  6.  Structure  with  clamped-clamped  boundary  conditions 
governing  the  transverse  vibrations  and  clamped-sliding  boundary 
conditions  governing  the  axial  vibrations. 

For  later  use,  the  following  dimensionless  quantities  are 


(4) 


For  the  microresonator  parameter  values  given  in  table  3,  the 
slenderness  ratio  r  =  y/EI /EA  l~l  is  of  the  order  of  10-3. 
Based  on  this  observation,  an  additional  assumption  that  U  = 
0(W2)  has  been  used,  where  Of)  is  the  Landau  symbol  that  is 
used  here  to  state  that  U  is  of  the  order  of  W2  [18].  It  is  also 
assumed  that  W  =  O(r)  and  e  =  O^r2). 

Introducing  the  dimensionless  quantities  given  by 
equation  (4)  into  equations  (l)  and  (2),  dropping  the  higher- 
order  terms,  and  simplifying,  the  governing  equation  of  a 
uniform  beam  with  constant  pA ,  El  and  EA  values  can  be 
obtained  as  shown  below. 


EI(x)  =  YJEhHk(x),  /i(jr)  =  £>*tf*(i), 

A=1  *= I 

where  the  subscript  k  represents  the  klh  section  of  the 
resonator  and  this  index  runs  from  1  to  3.  The  quantities 
pAk ,  Elk ,  EAk  and  hk  are  constant  over  the  section  of 
interest.  Furthermore,  in  equation  (8),  the  function  //*(i)  = 
[u(x  -  Xk~\)  -  u(x  —  Jc*)]  is  constructed  by  using  the  unit 
step  function  u(x). 

Following  the  earlier  work  on  structures  with  piezoelectric 
elements  (e.g.,  [21]),  noting  that  the  distributed  actuation  is 
present  in  the  first  section  of  the  microresonator,  the  distributed 
loading  is  determined  as  (figure  2) 

ft(xj)  =  p(t)&(x  -x\),  p(t)  =  pq  cos  Clt, 

-  -  .  r  .  .  (9) 

where  p0  is  the  input  forcing  amplitude,  Cl  is  the  input 
excitation  frequency  and  <$(Jc)  is  the  Dirac  delta  function. 
For  convenience  of  discussion,  the  following  quantities  are 
introduced: 

EAk 


Xk  =  Uk(xk), 


EA  = 


[M 


Kk  = 

( xk  -xk-\) 

N0  =  ~EAe 


n 

r  _ 

n 

EA]  = 

- 1 

7.» 

* 

_ i 

ea2  = 

E  «r' 

_i=m- fl 

(10) 


where  ()  means  an  average  value. 

After  using  the  boundary  conditions  given  in  (7),  the 
displacement  compatibility  conditions 

Uk(xk)  =  UkMxk),  Wk(xk)  =  W*+|  (**)  (11) 

and  making  use  of  equation  (9)  in  (6),  the  result  obtained  is  the 
following  equation  governing  the  transverse  vibration  field  in 
the  kih  section  [19]: 

32W*  aW*  ekd2Wk  i)4W*  .  „ 

+  ~  7*S*  +  1*  =  FkC0S°klk 

xe[xk.uxt),  k  =  1,2,3.  (12) 

In  equation  (12), 


ek{tk)  =  rlNk+rlqkpk  cos (fi*f*)  + 


x  e  [**_i 
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Aik  = 
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EA 
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Tax' 
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k  ^  m 
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(13) 


(14) 


(15) 
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Nk  = 


Pk  = 


No/2 
Elk  ' 

Pol2 

Eh’ 


EkW  =  Pkhk 


35(jc  —  jci) 
dx 


and 


V  Elk 


(16) 


Equations  ( 1 2)— <  1 6)  are  suitable  for  studies  of  structures 
with  axially  varying  properties  and  distributed  actuation.  In 
addition,  this  model  is  applicable  to  structures  with  ‘small’ 
slenderness  ratios  such  as  the  microresonators  considered  in 
this  work. 

In  sections  3.3,  4  and  5,  solutions  of  equation  (12)  are 
sought.  First,  the  possibility  for  static  buckling  is  explored  and 
this  corresponds  to  the  situation  when  the  resonator  structure 
is  released  from  the  Si  wafer.  Following  that,  free  and  forced 
oscillation  problems  around  the  post-buckling  state  are  then 
solved.  The  developments  of  sections  4  and  5  are  applicable 
to  any  resonator  with  a  non-flat  equilibrium  position. 


3.3.  Static  buckling  problem 

In  this  section,  the  static  buckling  problem  is  addressed  along 
the  lines  of  the  authors’  previous  work  [17].  It  is  assumed  that 
the  dimensionless  static  buckling  deflection  of  the  resonator 
can  be  expressed  as 

n 

X(x)  =  J2xk(x)Hk(x),  (17) 

*  =  1 

where  Xk  (•*)  is  the  transverse  displacement  of  the  k\ h  section  of 
the  beam.  For  the  static  buckling  analysis,  the  inertia,  damping 
and  external  actuation  terms  in  equation  (12)  are  dropped, 
and  the  static  axial  forces  (extensional)  N0  and  Nk  given  in 
equations  (10)  and  (16)  are  replaced  by  the  static  compressive 
forces  —  Po  and  —  Po,k,  respectively.  After  carrying  these 
substitutions  through,  the  nonlinear  equation  governing  the 
equilibrium  of  the  beam  subjected  to  a  static  axial  force  can 
be  written  as 


*=  1,2,. (18) 

The  associated  boundary  conditions  follow  the  form  of  (7). 

The  approximate  solution  of  the  post-buckling  system 
described  by  equation  (18)  can  be  obtained  as  given  below, 
in  terms  of  a  critical  buckling  mode  shape  with  the  amplitude 
b ,  called  the  buckling  factor. 

X/t(x)  =  bif/kix).  (19) 


Hi,i)  Dynamic  displacement  M*J)  :  Total  displacemant 


Figure  7.  Buckled  beam  configuration. 


by  making  use  of  the  following  boundary  conditions  and 
compatibility  conditions: 


=  o,— - =0 

At 

fn(x)  =  0,— - - =0 

dx 

'I'k-iM  =  \l>k(x), 


at  x  =  xo, 

at  *=*„. 

dV^-iOO 


dx 


djc 


d2^»-i(x) 

d V»(x) 

Eh- 1 

1  dx2 

dV*-i(x) 

Eh- 1 

1  dr3 

- 

at  x  =  xk- 

(22) 


(23) 


Here,  it  is  important  to  note  that  in  making  the  approximation 
(19),  the  boundary  condition  for  axial  motions  is  held  fixed 
at  x  =  /  instead  of  the  sliding  boundary  condition  shown  in 
figure  6,  previously. 


4.  Free  linear  oscillations  about  the  post-buckled 
position 


Free  oscillations  of  the  undamped  resonator  are  considered  in 
this  section.  The  governing  equation  of  motion  in  the  presence 
of  a  compressive  axial  load  is  of  the  form 


d2wk 

dx2 


34l Vk 

+ - f=0,  *=1,2 . n.  (24) 

dx4 

where  wk  is  the  overall  transverse  displacement  in  the  free- 
vibration  case.  The  corresponding  boundary  and  compatibility 
conditions  are  similar  to  those  given  by  equations  (22)  and  (23). 
Following  the  authors’  previous  work  ([17,  18]),  the  solution 
of  this  boundary-value  problem  can  be  written  as  the  sum  of  a 
static  component  and  a  dynamic  component  (see  figure  7). 


Wk(x,  tk)  =  b\!/k{x)  +  vk(x ,  /*),  (25) 


This  is  given  by 

yi  _  _ (Po  ~  Pc) _  (20) 

EAk  ET-I  [A*»  fxi,  (d'Ai/Ar)2^]’ 

where  \)/k  ( x )  is  the  *th  component  of  the  static  critical  buckling 
mode  shape  yj/(x)  and  it  has  the  form 

\//k(x)  =aLk  +  al  kx  +  a3.*  cos(fr*)  +  aA  k  sin(£**).  (21) 

In  equation  (21),  =  Pc  k  =  Pcl2/Elk  and  the  constants  a,  * 

for  i  =  1 , 2,  3,  4  and  k  =  1 , 2, . . . ,  n  can  be  determined 


where  vk(x ,  tk )  is  the  *th  component  of  the  dynamic  deflection. 

Next,  the  natural  frequencies  cb  and  mode  shapes 
associated  with  free  oscillations  around  the  post-buckled 
position  b\//(x)  are  examined.  To  find  the  natural  frequencies 
and  mode  shapes,  one  first  substitutes  equation  (25)  into 
equation  (24),  assumes  that  |u*(jc,  f*)|  <K  16^*001,  and  uses 
separation  of  variables  in  the  following  form: 

vk(xttk)  =  Mx)t^\  (26) 

where  j2  =  —1,  and  <pk(x)  and  cok  =  pAk/ Elkl2a)  are 
the  mode  shape  and  the  natural  frequency  associated  with 
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section  k ,  respectively.  Here,  it  is  necessary  to  emphasize  the 
relationship  between  a>*  and  u)\  that  is,  a>*f*  =  cot.  Making 
use  of  equations  (20),  (24),  (25)  and  (26),  the  following  is 
obtained: 


2  d24>k 

-d)2k<t>k  +  Pck~^2  + 


dV* 

d*4 


2 b2  d  V* 
rl  ^ 


The  general  solution  of  (27)  is  given  by 


(27) 


<PkW  =  C,.*  sin(Xi.*jr)  +  C2.k  cos(Xux)  +  C3*  sinh(A2,*jr) 


+  C41  cosh(A2,**)  +  C5jt[a3.*  cos((kx)  +  a4.t  sin«*jt)], 


(28) 

where 

_  (29) 

A2.*  =  >/K-^.‘  +  V/pc2*+4"r) 
and  the  Cr*,  r  =  1 , 2, . . . ,  5,  are  constants  that  satisfy 

r  A,*(^i.,Ci.i  +  B2jCi  i  +  Bi  jCij  +  BijCt  i  +  DkjCyi) 

i=i 

=  0,  (30) 

where  the  £r>l  and  D*<f  (r  =  1, 2,  3,  4  and  k  =  1, 2, . . . ,  n) 
are  constants,  the  expressions  for  which  are  detailed  in 
[17,  19]  along  with  the  derivations  for  equations  (24)-(30). 

Equation  (30)  along  with  the  associated  boundary 
conditions  and  compatibility  conditions  forms  an  eigenvalue 
problem  for  the  eigenmodes  Cr>,  and  the  natural  frequencies 
(Dk.  After  determining  these  values  for  Cr>(  ,  the  expression  for 
the  free  vibration  mode  shapes  of  the  post-buckled  beam  with 
axially  varying  properties  can  be  determined  as 

n 

<t>(x  )  =  £&(*)//*(*).  (31) 

k=\ 


5.  Forced  oscillations 


Forced  oscillations  of  a  resonator  modeled  as  a  beam  with 
axially  varying  properties  are  considered  in  this  section.  The 
equation  of  motion  is  given  by  equations  (12)  and  (13), 
the  boundary  conditions  are  given  by  equations  (7),  and 
the  displacement  compatibility  conditions  are  similar  to  those 
given  by  (23).  For  solving  (12),  it  is  assumed  that 


W*(*.  tk)  =  br/rk(x)  +  rjk(x,  tk),  (32) 


where  77* (jc,  tk)  represents  the  dynamic  deflection  and  b\j/k(x) 
is  the  static  post-buckling  deflection  given  by  ( 1 9).  Making  use 
of  the  development  given  in  section  3.3  and  (32),  equation  ( 1 8) 
cart  be  rewritten  as 


„  in 

i)k  +  Pc.kT)k  +*)k  ~  ~ T 


hikWi,  vii)i 


=  [Fk+qkpkWi  +  »J*)]cos(Qj/*)  -  2fikrik 
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r*  L.=i  J 

■?|e 

r*  Li-i 


rh  +  ~2 


Y2  v'j)i 


rh 


,  A,»(b,'.  nt)t 


k  =  1,2 . n.  (33) 


From  here  onwards,  the  following  notation  and  definitions  are 
used  for  convenience:  ()  =  d/d  tk,  ()'  =  d/dx,  (a,b)j  = 
f*‘  j  a(x)b(x)  dx  and  (a,  b)  =  f**  a(x)b(x)  dx.  Again,  the 
boundary  and  compatibility  conditions  are  similar  to  those 
given  by  (22)  and  (23).  Next,  assuming  weak  damping  and 
forcing  terms,  the  method  of  multiple  scales  [22]  is  used 
to  seek  for  an  approximate  solution  of  (33).  The  basic 
approach  follows  that  presented  in  an  earlier  work  [16],  but 
here,  the  application  is  directed  toward  a  composite  microscale 
structure  with  varying  axial  properties.  Different  time  scales 
are  introduced  by  using  a  small,  dimensionless  book-keeping 
parameter  e  as  follows: 

Tok  =  tk,  T\k  =  etk ,  72*  =  e2tk,  ....  (34) 


The  time  derivatives  take  the  form 

d/dtk  =  Dok  +  £ D\k  +  e2 Dik  +  •  •  - ,  Dik  =  d/dTjk, 

(35) 

where 


To  balance  the  effect  of  the  nonlinearity,  damping  and 
excitation,  the  following  scaling  is  used: 

Pk  =  e3Pk,  A*  =  £2P-k-  (37) 


The  approximate  solution  is  then  expanded  as 

Vk(x,  tk)  =  erjik(x,  T0k,  T\k,  Tik)  +  £2*l2k(x,  Tok ,  T\k,  Tzk) 

+  £3ri3k(x,  7o*.  T\k,  Tu)  +  •  •  •  (38) 


Introducing  (34)-(38)  into  (33),  rewriting  the  excitation  term 
in  polar  form,  and  collecting  terms  of  the  same  power  of  e,  the 
following  hierarchy  of  equations  is  obtained: 

0(8): 

C(7Uk)  =  Dlkrik  +  Pc.kq"k  +  tj\\ 


0(e2): 
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where  c.c.  stands  for  the  complex  conjugate  of  the  preceding 
quantity  and 

HS(x  -xm) 

Pk  —  Pkbk - - - +  bqkpk\//k .  (42) 

dX 
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Figure  8.  Results  of  nonlinear  analysis:  ( a )  first  natural  frequency  versus  b ,  ( b )  predicted  forced  spatial  response  when  b  =  4.65  x  10  4  and 
( c )  frequency-response  curve  for  b  =  4.65  x  10-4,  p  =  4.65  x  10-5  and  £  =  8.23  x  10-4.  The  solid  and  dashed  lines  represent  stable  and 
unstable  branches,  respectively. 


The  corresponding  boundary  and  compatibility  conditions  at 
all  orders  are  similar  to  those  given  by  equations  (22)  and 
(23).  The  excitation  frequency  is  assumed  to  be  close  to  the 
first  natural  frequency;  that  is, 

Q  =  &i  +  e2(7,  or  Qk  =  o)\k  +  e2Ok-  (43) 

The  solution  of  (39)  can  be  written  as 

riik  =  0uW[^.(r,*.r2*)e^»r“+^1(ru,r2*)e-^‘r“]. 

(44) 

where  A  i  is  the  complex  conjugate  of  A  j ,  which  can  be  written 
in  polar  form  as  A\{T\k ,  7^)  =  On  substituting  (44) 

into  (40),  the  right-hand  side  of  (40)  can  be  expressed  as  the 
summation  of  terms  that  produce  non-secular  terms  and  secular 
terms  (ST), 

C(V2k)  =  fc4'1W[<4?ej2“'"7"  +  A\A\  +c.c.]  +  ST,  (45) 


where 


».w- 4  £>(*;•«!>» 

r‘  Li-i  J 
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ST=  — j2«u^*(D|*A1ej"-‘7'“  -  D,*A,e- 


(46) 


Setting  the  sum  of  the  secular  terms  to  zero  one  obtains 
A  j  =  A\{Jik ),  and  solving  (45)  by  assuming  the  approximate 
solution  for  the  spatial  term  of  rj2k(x,tk)  to  be  a  weighted 
summation  of  the  free-oscillation  mode  shapes  given  by  (28), 
one  has 

tiu  =M,A,<l>u(x)  +  M;<D2*U)ej2“,*7'“  +c.c„  (47) 


where  4>i*(x)  and  4>2a(x)  are  defined  as 

oo 

4>u(x)  =  ^Trfrk, 

r=l 
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Tr  = 
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(48) 
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(49) 

EL,  EL,  W.  4-  M-  trth] 

l2[6>2-(  2&,)2]<0ri0,> 

Next,  to  solve  (41),  substituting  (43),  (44)  and  (47),  one  arrives 
at 

£(«J3»)  =  *?*,*(*)  eJ*““r"  +hu(x,  ra)e^‘r“  +c.c., 


where 
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2b2 


solution  of  the  corresponding  homogeneous  adjoint  problem 
from  x  =  0  to  1,  which,  from  (26)  and  (31),  is  <f>rkW  c±i&rt 
for  the  kt h  section  and  the  rth  mode.  Multiplying  the 
right-hand  side  of  (50)  by  <f>\k(x)  spatially  integrating 

the  result  from  x  =  jc*_i  to  jc  =  adding  all  sections 

together,  and  setting  the  sum  of  the  secular  terms  to 
zero,  the  following  complex-valued  modulation  equation  is 
obtained: 


[y  dukA  i + D2kA\)(v\k(p\k,  <p\k)k 

k=) 


y  (gk,  </>u)* 


(50) 


(54) 

Further,  substituting  (4)  and  (36)  into  (54),  and  letting 
y  =  (7  7*2  —  P,  it  is  found  that  the  approximate  solution  and 
frequency-response  equation  are  respectively  given  by 
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The  solvability  condition  [22]  requires  the  right-hand  side 
of  (50)  and  its  boundary  conditions  be  orthogonal  to  every 


where 
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Figure  9.  Comparisons  between  experimental  measurements  and  the  predicted  frequency  responses:  ( a )  thickness  is  90  nm, 

£  =  8.23  x  10~4,  b  =  4.65  x  10-4  and  ann  =  -2.73  x  10";  ( b )  thickness  is  93.77  nm,  £  =  8.23  x  10"\  b  =  4.65  x  10'4  and 

&nn  =  -2.84  x  1010;  (c)  thickness  is  140  nm,  £  =  8.23  x  10~4,  b  =  4.85  x  10~4  and  ann  =  -2.77  x  10n;  and  ( d )  thickness  is  160  nm, 

£  =  8.23  x  IQ'4,  b  =  4.86  x  10~4  and  ann  =  4.36  x  1010. 


6.  Results  and  discussion 

In  figure  8(a),  the  variation  of  the  first  natural  frequency  of  a 
200  /zm  PZT  resonator  versus  the  buckling  factor  b  is  shown 
for  free  oscillations  about  the  fifth  static  buckling  mode.  It 
is  important  to  note  that  for  the  same  frequency  value,  the 
corresponding  value  of  b  depends  on  the  resonator’s  stiffness 
and  geometric  properties  as  well  as  the  normalization  used 
to  define  the  static  critical  buckling  mode  shape.  (Here, 
the  normalization  used  for  the  static  critical  buckling  mode 
shape  \//(x)  is  /J  ^2(jr)djc  =  1.)  When  b  equals  4.65  x 
10-4,  the  first  natural  frequency  of  the  resonator  is  determined 
to  be  3 1 3  kHz,  which  is  close  to  the  experimental  result  given  in 
the  second  section  of  this  paper.  The  corresponding  deflection, 
b[\J/  +  \ba2$> i]  +  a<p i,  determined  from  equations  (32)  and 
(44),  is  shown  in  figure  8 (b).  In  this  figure,  the  amplitude  a 
of  mode  shape  <p  is  scaled  to  be  10-3,  so  that  the  predicted 
dynamic  displacement  is  of  the  same  order  as  that  seen  in 
experiments.  When  the  excitation  frequency  is  close  to  the 
first  natural  frequency  of  the  resonator,  the  predicted  forced 


spatial  response  is  similar  to  the  experimentally  observed 
spatial  pattern  (see  figure  4(a)).  The  frequency-response 
curve  obtained  from  (55)  is  shown  in  figure  8(c).  It  resembles 
that  of  a  Duffing  oscillator  with  a  hardening  spring  and  this 
is  qualitatively  consistent  with  the  experimental  results  of 
figure  3(b).  Guided  by  related  work  [23],  the  damping 
coefficient  £  and  the  excitation  force  amplitude  p  are  chosen 
to  generate  the  predictions. 

The  agreement  between  the  analytical  predictions  and 
experimental  data  suggests  that  the  hypothesis  of  the  non¬ 
flat  equilibrium  position  of  the  resonator  caused  by  buckling 
can  be  a  valid  one.  The  present  work  can  be  used  as  a 
basis  to  study  buckling  and  oscillations  of  microresonators 
with  axially  stepwise  varying  properties.  For  different  values 
of  the  top  electrode  thickness,  the  numerical  values  of  ann 
and  C  are  found  to  vary.  The  associated  comparisons 
between  the  experimentally  observed  frequency  response  and 
the  predicted  frequency  response  obtained  from  equation  (55) 
are  given  in  figure  9.  The  dotted  lines  correspond  to  the 
experimental  measurements  and  the  solid  lines  correspond 
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to  model  predictions.  It  is  seen  that  the  top  electrode 
thickness  affects  the  frequency-response  curve  dramatically 
by  changing  it  from  hardening  to  softening  in  some  cases. 

Acknowledgments 

The  authors  gratefully  acknowledge  the  support  received 
for  this  work  through  DARPA  contract  no  F3060202C00 1 6, 
AFOSR  grant  no  F49620-03-10181,  and  ARO  grant  no 
W91 1NF05 10076.  Material  property  data  were  obtained  from 
researchers  at  the  Maryland  MEMS  Laboratory,  University 
of  Maryland,  College  Park,  Maryland  (MD),  and  the  Army 
Research  Laboratory,  Adelphi,  MD. 

References 

[  1  ]  Lin  L,  Nguyen  C  T-C,  Howe  R  T  and  Pisano  A  P  1 992 
Microelectromechanical  filters  for  signal  processing 
Technical  Digest,  IEEE  Micro  Electromechanical  Systems 
Workshop  (Travemunde,  Germany,  4-7  February  1992) 
pp  226-31 


[2]  Younis  M  I  and  Nayfeh  A  H  2003  A  study  of  the  nonlinear 

response  of  a  resonant  microbeam  to  an  electric  actuation 
Nonlinear  Dyn.  31  91-1 17 

[3]  Balachandran  B  and  Li  H  2003  Nonlinear  phenomena  in 

microelectromechanical  resonators  Proc.  IUTAM  Symp.  on 
Chaotic  Dynamics  and  Control  of  Systems  and  Processes  in 
Mechanics  (Rome,  Italy,  8-13  June)  pp  97-106 

[4]  Piekarski  B,  DeVoe  D,  Dubey  M,  Kaul  R,  Conrad  J  and  Zeto  R 

2001  Surface  micromachined  piezoelectric  resonant  beam 
filters  Sensors  Actuators  A  91  3 1 3-20 

[5]  DeVoe  D  L  2001  Piezoelectric  thin  film  micromechanical 

beam  resonators  Sensors  Actuators  A  88  263-72 

[6]  Currano  L  2002  Experimental  and  finite  element  analysis  of 

piezoelectrically  driven  MEMS  actuators  MS  Thesis 
Department  of  Mechanical  Engineering,  University  of 
Maryland,  College  Park 

[7]  Anderson  T  J,  Nayfeh  A  H  and  Balachandran  B  1996  Coupling 

between  high-frequency  modes  and  a  low  frequency  mode: 
theory  and  experiment  Nonlinear  Dyn.  11  17-36 

[8]  Li  H  and  Balachandran  B  2002  Buckling  induced  nonlinear 

phenomenon  in  a  micro-electromechancial  resonator  Proc. 
ASME  Int.  Mechanical  Engineering  Congress  and 
Exposition  (Orlando,  FL,  2002)  Paper  No 
IMECE2002-3901 0 

[9]  Preidikman  S,  Li  H  and  Balachandran  B  2003  Forced 

oscillations  of  microelectromechanical  resonators  Proc . 
Symp.  on  Nonlinear  Dynamics  and  Stochastic  Mechanics, 


366 


Nonlinear  free  and  forced  oscillations  of  piezoelectric  microresonators 


2003  ASME  Int.  Mechanics  Engineering  Congress  and 
Exposition  (Washington,  DC,  15-21  November)  Paper  No 
IMECE2003-44552  [17] 

[10]  Li  H  and  Balachandran  B  2003  Nonlinear  oscillations  of 

micromechanical  oscillators  Proc.  2003  ASME  Int.  Design 
Engineering  Technical  Conf  ( Chicago )  Paper  No  [18] 

DETC2003/VIB-48520 

[11]  Chiao  M  and  Lin  L  2000  Self-buckling  of  micromachined  [19] 

beams  under  resistive  heating  J.  Microelectromech. 

Syst.  9  146-51 

[12]  Lindberg  U,  Soderkvist  J,  Lammerink  T  and  Elwenspoek  M 

1998  Quasi-buckling  of  micromachined  beams  [20] 

J.  Micromech.  Microeng.  3  183-6 

[13]  Fang  W  and  Wickert  J  A  1994  Post-buckling  of  [21] 

micromachined  beams  J.  Micromech.  Microeng.  4  1 16-22 

[14]  Saif  M  T  A  2000  On  a  tunable  bistable  MEMS — theory  and 

experiment  J.  Microelectromech.  Syst.  9  157-70  [22] 

[15]  Nayfeh  A  H  and  Lacarbonara  W  1998  On  the  discretization  of 

spatially  continuous  systems  with  quadratic  and  cubic  [23] 

nonlinearities  JSME  Int.  J.  C  41  5 1 0-3 1 

[16]  Lacarbonara  W  and  Nayfeh  A  H  1998  Experimental 

validations  of  reduction  methods  for  nonlinear  vibrations  of 


distributed-parameter  systems:  analysis  of  a  buckled  beam 
Nonlinear  Dyn.  17  95-1 17 

Li  H  and  Balachandran  B  2005  Buckling  and  free  oscillations 
of  composite  microresonators  J.  Microelectromech.  Syst. 
at  press 

Nayfeh  A  H  and  Mook  D  T  1979  Nonlinear  Oscillations  (New 
York:  Wiley)  pp  447-53 

Li  H  2006  Oscillations  of  microscale  composite  structures 
with  applications  to  microresonators  PhD  Dissertation 
Department  of  Mechanical  Engineering,  University  of 
Maryland,  College  Park 

Talpaert  Y  R  2002  Tensor  Analysis  and  Continuum  Mechanics 
(Dordrecht:  Kluwer)  pp  199-209 

Crawley  E  F  and  Lazarus  KB  1991  Induced  strain 
actuation  of  isotropic  and  anisotropic  plates  AIAA  J. 

29  944-51 

Nayfeh  AH  1981  Introduction  to  Perturbation  Techniques 
(New  York:  Wiley)  pp  388^*71 

Dick  A  J,  Balachandran  B,  DeVoe  D  L  and  Mote  C  D  Jr  2005 
Parametric  identification  of  piezoelectric  micro-scale 
resonators  Proc.  ENOC-2005  ( Eindhoven ,  Netherlands, 
7-12  August  2005) 


367 


Institute  of  Physics  Publishing 


Journal  of  Micromechanics  and  Microengineering 


J.  Micromech.  Microeng.  16(2006)512-525  doi:  1 0. 1 088/0960- 131 7/1 6/3AXJ6 

A  semi-analytical  tool  based  on  geometric 
nonlinearities  for  microresonator  design 

Sergio  Preidikman  and  Balakumar  Balachandran 

Department  of  Mechanical  Engineering,  University  of  Maryland,  College  Park, 

MD  20742-3035,  USA 

Received  22  September  2005,  in  final  form  23  December  2005 

Published  30  January  2006 

Online  at  stacks.iop.org/JMM/16/512 

Abstract 

In  this  paper,  a  computational  mechanics  model  specifically  tailored  for 
composite  microresonators  with  piezoelectric  actuation  and  piezoelectric 
sensing  is  developed  and  used  as  a  design  tool  for  these  microresonators. 

The  developed  model  accounts  for  the  structural  properties  and  the 
electromechanical  coupling  effect  through  finite-element  analysis.  It  is 
assumed  that  the  deflection  is  large  and  that  the  geometric  nonlinearity  must 
be  included.  The  dynamic  admittance  model  is  derived  by  combining  the 
linear  piezoelectric  constitutive  equations  with  the  modal  transfer  function 
of  the  multi-layered  microresonator  structure.  The  resonator  receptance 
matrix  is  constructed  through  modal  summation  by  considering  a  limited 
number  of  dominant  modes.  The  electromechanical  coupling  determination 
at  the  input  and  output  ports  makes  use  of  converse  and  direct  piezoelectric 
effects.  In  the  development  of  the  finite-element  models,  the  boundary 
conditions,  the  shapes  of  electrodes  and  distributed  parameters  such  as 
varying  elastic  modulus  across  the  length  of  the  structure  have  been  taken 
into  account.  The  developed  semi-analytical  tool  can  be  used  to  carry  out 
parametric  studies  with  respect  to  the  following:  (i)  the  resonator  beam 
thickness  and  length,  (ii)  the  influence  of  constant  axial  forces  on  the 
transverse  vibrations  of  clamped-clamped  microresonators,  (iii)  the 
geometry  of  the  drive  and  sense  electrodes  and  (iv)  imperfect  boundary 
conditions  due  to  mask  imperfections  and  fabrication  procedure.  The 
semi-analytical  development  has  been  validated  by  comparing  model 
predictions  with  prior  results  available  in  the  literature  for  clamped-clamped 
resonators  and  experimental  measurements.  A  detailed  discussion  of 
modeling  considerations  is  also  presented. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Piezoelectrically  actuated  microscale  resonators  are  attractive 
for  communication  and  signal-processing  applications  [1]. 
Two  types  of  resonators  are  considered  here.  One  type  of 
resonators  will  be  referred  to  as  the  AlGaAs  resonator,  and 
the  other  type  of  resonators  will  be  referred  to  as  the  PZT 
resonator.  Both  types  of  resonators  are  composite  structures, 
and  the  PZT  resonators  have  asymmetric  cross-sections,  as 
discussed  in  previous  work  [2,  3].  The  considered  resonators 
are  based  on  the  piezoelectric  effect,  as  shown  in  figure  1. 
The  elastic  substrate  is  a  SiC>2  layer,  on  the  top  of  which  a 
platinum  electrode  layer  is  deposited  throughout  the  length 


of  the  structure.  A  thin  layer  of  sol-gel  piezoelectric  film  is 
located  on  the  top  of  this  electrode  layer.  To  complete  the 
structure,  another  platinum  layer  is  deposited  on  the  top  of 
this  piezoelectric  film  and  this  layer  extends  over  one  quarter 
of  the  length  from  each  anchor.  The  mid-section  of  the 
resonator  structure  is  free  from  this  platinum  electrode  layer 
[1].  Due  to  the  asymmetry  of  the  cross-section,  the  position 
of  the  piezoelectric  layer  is  offset  from  the  neutral  axis,  and  in 
addition,  (tensional)  residual  stress  may  also  be  introduced  in 
each  layer  during  the  fabrication  process.  The  effect  of  these 
stresses  is  explored  in  this  work. 

As  pointed  out  in  the  authors’  recent  work,  the 
resonators  also  exhibit  nonlinear  characteristics  [4-6].  These 
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Figure  1.  ( a )  SEM  of  a  PZT  resonator  (courtesy,  Maryland  MEMs 
Laboratory)  and  ( b )  a  schematic  showing  the  details  [2]. 


Suspended  Beam  Substrate  Piezoelectric  Film 


Figure  2.  A  clamped-clamped  composite  resonator  subjected  to  a 
constant  axial  load  P. 

characteristics  include  Duffing  oscillator-like  response  during 
resonance  excitations,  temporal  harmonics  in  the  response 
and  spatial  patterns  during  forced  oscillations  that  cannot  be 
explained  by  a  conventional  linear  analysis.  The  lengths  of 
the  resonators  considered  in  previous  studies  typically  range 
from  100  /xm  to  400  /xm,  and  the  thickness  of  each  platinum 
electrode  is  in  the  range  of  90  nm  to  180  nm  [2-6].  In  some 
typical  uses  of  this  resonator,  the  structure  is  driven  close  to  its 
first  resonance  frequency  with  the  input  at  the  drive  electrode 
having  a  dc  bias  in  addition  to  the  harmonic  component. 

In  this  work,  the  authors  discuss  a  semi-analytical  finite- 
element  (e.g.  [7,  8])  based  formulation,  in  which  transverse- 
free  vibrations  of  composite  and  axially  stepwise  varying 
properties  of  microresonators  subjected  to  constant  axial  loads 
are  considered  (see  figure  2  for  a  clamped-clamped  case). 
It  is  shown  that  the  consideration  of  axial  loads  is  important 
to  predict  the  natural  frequencies  of  the  resonators  observed  in 
the  experiments. 

To  take  into  account  the  effect  of  the  axial  load,  a 
geometrically  nonlinear  analysis  needs  to  be  performed.  In 
such  a  case,  the  equations  of  motion  for  the  microresonator 
can  be  written  in  the  following  form: 

Md(0  +  Cd(0  +  Kd(r)  =  f(0. 


where  K  =  Kf  +  Kg  is  the  global  stiffness  matrix,  is  the 
elastic  stiffness  matrix  and  Kg  is  the  geometrical  stiffness 
matrix  that  is  obtained  from  the  nonlinear  component  of 
the  strain-displacement  relation.  This  matrix  has  not  been 
considered  in  the  previous  studies  of  microelectromechanical 
systems  (MEMs). 

In  the  last  several  years,  there  has  been  a  growing  need  for 
accurate  modeling  and  simulation  of  microelectromechanical 
devices  and  systems  that  employ  piezoelectric  materials.  This 
comes  from  the  need  to  reduce  design  iterations  and  speed 
up  the  product  development,  and  also  to  ensure  reliability 
of  the  final  product.  Finite-element  analysis  (FEA)  plays  an 
important  role  in  the  simulation  of  MEMs  devices,  and  this 
analysis  generally  covers  multiple  domains  for  a  single  device, 
such  as  structural,  thermal,  electrostatic,  electromagnetic 
and  fluid  domains,  von  Preissig  and  Kim  [9]  examined 
techniques  for  modeling  thin-piezoelectric  MEMs  devices 
by  using  existing  finite-element  packages.  In  this  work, 
piezoelectrically  actuated  bending  is  examined.  The  authors 
point  out  that,  while  it  may  seem  that  the  sheet-like  nature 
of  structures  in  piezoelectric  MEMs  would  make  them  good 
candidates  for  conducting  FEA  with  plate  elements,  solid  or 
‘brick’  elements  can  work  remarkably  well.  Finite-element 
model  (FEM)  errors  associated  with  the  discretization  of  the 
model  have  also  been  analyzed.  An  important  issue  to  note 
is  that  meshing  a  thin  sheet  into  low-aspect-ratio  elements 
requires  a  prohibitively  large  number  of  elements,  while  too 
low  a  mesh  density  might  result  in  severe  discretization  and 
element-shape  errors.  A  four-node,  isoparametric,  linear 
piezoelectric,  plane-strain  element  from  the  ANSYS  library 
has  been  used  in  this  work.  Wang  and  Ostergaard  [10]  used 
finite  elements  to  develop  a  coupled  simulation  method  for 
piezoelectric  transducers  with  an  attached  electric  circuit. 
In  this  work,  the  weak  form  of  the  laws  of  conservation 
of  momentum  and  electric  charge  for  a  linear  piezoelectric 
medium  are  discretized  by  using  the  FEM.  Their  method 
has  been  implemented  in  the  ANSYS  software.  Chen  et  al 
[11]  presented  a  two-dimensional  analytical  model  of  a  spiral¬ 
shaped  PZT  ceramic  actuator.  They  used  the  FEA  to  validate 
the  results  obtained  from  the  analytical  model.  In  this  work, 
the  commercially  available  software  packages  PATRAN  and 
ABAQUS  are  used.  PATRAN  is  utilized  as  a  pre  processor 
and  ABAQUS  is  used  as  a  post-processor  to  perform  the  linear 
elastic,  piezoelectric  analysis.  In  order  to  capture  bending 
effects  accurately,  the  authors  used  eight-node,  isoparametric, 
plane  strain,  linear  elastic,  piezoelectric  elements.  For 
achieving  convergence,  while  keeping  the  length-to-width 
ratio  of  the  elements  reasonable,  they  used  meshes  with  at 
least  ten  elements  across  the  spiral  thickness.  A  typical  finite- 
element  mesh  of  a  two-tum  spiral  actuator  has  more  than 
3000  elements. 

In  this  paper,  the  coupling  of  electrical  and  mechanical 
fields  that  is  intrinsic  to  a  piezoelectric  material  is 
accomplished  in  a  non-traditional  approach.  From  a  filter 
design  standpoint,  the  admittance  function  relating  input 
voltage  to  output  current  l2{(v)/ V\{co)  is  an  important 
frequency-response  function  to  be  determined  (see  figure  3). 
This  admittance  function  may  be  obtained  by  relating  the 
mechanical  transfer  function  of  the  microresonator  structure 
to  the  corresponding  electrical  input  and  output  through  the 
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Figure  3.  Input  and  output  quantities  of  interest  for  a  clamped- 
clamped  piezoelectric  microresonator. 

piezoelectric  constitutive  equations.  In  the  second  section 
of  this  work,  this  admittance  function  is  determined  for 
a  composite  microresonator  with  axially  stepwise  varying 
properties.  This  microresonator  is  also  subjected  to  constant 
axial  loads.  The  obtained  admittance  function  has  the  form 

K2,(o>)  =  -£7^  =  C*(o>)Hm(o>)F\ 

where  H5A(o>)  is  a  ‘block’  of  the  complex  frequency  response 
matrix  for  the  considered  microresonator  (e.g.  [7]);  C s(co)  is 
a  matrix  that  depends  on  the  material  properties  of  the  PZT 
material,  the  geometry  of  the  sensor  and  the  finite-element 
discretization  of  the  microresonator,  and  F*  is  a  vector  that 
depends  on  the  material  properties  of  the  PZT,  the  geometry  of 
the  actuator,  the  geometry  of  the  suspended  beam  substrate  and 
the  finite-element  discretization  of  the  microresonator.  In  the 
second  section  and  the  appendices,  the  expressions  for  these 
matrices/vectors  are  derived.  Comparisons  with  experimental 
results  and  results  available  in  the  literature  are  presented  in  the 
third  section.  It  is  believed  that  the  numerical  and  analytical 
efforts  presented  in  this  work  can  be  used  as  a  basis  to  develop 
design  tools  for  piezoelectric  microresonators. 


where  r,  s  =  1,2,...,  n& Df,  fldof  is  the  number  of  degrees  of 
freedom  (dof)  of  the  finite-element  model,  a>i  is  the  /th  natural 
frequency  of  the  undamped  system  and  £/  is  the  /th  damping 
ratio  of  the  microresonator,  which  is  given  by 

</  =  r—  (or  +  /*«>/)•  (5) 

2  o)i 

In  (4),  (pir  and  <pis  are,  respectively,  the  r-component  and  the 
5-component  of  the  /th  mass  normalized  mode  shape  of  the 
microresonator  and  these  modes  satisfy 

Y’mM'A,  =  &mn.  (6) 

Further,  in  (4),  nmod  is  the  number  of  dominant  modes.  The 
natural  frequencies  and  mode  shapes  in  (4)  are  determined 
from 


K3>  =  M^A,  (7) 

where 

$  =  [y’i ¥>2  •••¥’*,**].  A  =  [<Sm„aj2].  (8) 

The  damping  ratio  in  (4)  is,  in  general,  determined 
experimentally  from  the  measured  quality  factor,  Qi ,  which  is 
given  by 


QCu  ttcl 


where  the  cutoff  frequencies  are  given  by  [  1 3] 


(9) 


2.  Semi-analytical  development 


2.7.  The  equations  of  motion  in  the  frequency  domain 


The  equations  of  motion  of  the  microresonator  in  the  frequency 
domain  are  of  the  form 

{— a>2M  +  ja>C  +  K}d(cu)  =  f(a;),  (1) 

where  j  =  V^T,  K  =  +  Kc  is  the  global  stiffness  matrix, 

M,  Kf  and  Kg  are  the  global  mass  matrix,  elastic  stiffness 
matrix  and  geometrical  stiffness  matrix,  respectively,  and  d(co) 
and  f(a>)  are  the  Fourier  transforms  of  the  nodal  displacements 
d(f)  and  the  nodal  forces  f(t),  respectively.  The  development 
of  the  time-domain  model  with  the  geometric  nonlinearity  is 
detailed  in  appendix  A.  From  ( 1 ),  one  obtains 

d(a>)  =  H(co)f((v),  (2) 


where 

H(a>)  =  {— ta2M  +  j<wC  +  K}~ 1  (3) 


is  the  complex  frequency  response  matrix  (the  mechanical 
force-displacement  transfer  function)  of  the  beam 
microresonator  [7],  Since  proportional  damping  is  assumed, 
the  complex  frequency  response  matrix  can  be  constructed 
from  the  modal  summation  as  [  12] 


H(o>)  =  [HrAu)]  = 


"mod 

E 


_/=i 


_ 4>lr<Pls _ 

0>?[1  -  (w/lD,)2  +j2C;(a>/<t>,)] 


(4) 


ncl  =  y]l  -  -  2<,/l  - 1,2, 

(10a) 

=  y/l  -  2 c,2  +  2C,/l  -  C2. 

(106) 

After  substituting  (10)  into  (9)  and  using  a  Taylor’s  series 
expansion  about  Qi  =  0,  the  result  obtained  is 

Qi  =  +  0(tf)-  (11) 

Hence,  for  very  lightly  damped  microresonators  (high-0 
microresonators),  that  is,  0  <  £  <£  1 ,  the  quality  factor  can  be 
approximated  as 


(12) 

Now,  in  terms  of  the  quality  factors,  it  is  possible  to  rewrite  the 
dof-to-dof  mechanical  complex  frequency  response  function 
of  the  microresonator  structure,  Hrs((o),  as 

^  d>,rs,, 

H„(o>)  =  Y  - .  (13) 

“ '  “>?U  -  (“/"/)  +j/£M"M)l 

At  this  point,  one  has  the  elements  of  the  complex  frequency 
response  matrix,  which  describes  the  resulting  harmonic 
displacement  of  the  rth  dof  caused  by  a  unit  harmonic  force 
applied  at  the  5th  dof. 
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Figure  4.  A  schematic  view  of  the  actuator  port. 


2.2.  Admittance  model 


In  this  section,  a  semi-analytical  model  describing  the 
admittance  function  for  a  clamped-clamped  piezoelectric 
microresonator  is  developed.  The  semi-analytical  admittance 
model  may  be  obtained  by  relating  the  (numerical)  mechanical 
transfer  function  of  the  microresonator  structure  to  the 
corresponding  electrical  input  and  output  through  the 
(analytical)  piezoelectric  constitutive  equations.  The  overall 
frequency-response  function  describing  the  admittance  Y2]  (<w) 
between  the  driven  input  (port  1)  and  the  sense  output  (port  2) 
of  the  beam  resonator  shown  in  figure  3  is  defined  as 


Y2\(co)  = 


h(co) 


(14) 


V\(co)' 

where  Vi(o>)  is  the  voltage  applied  to  the  input  port  and  I2(a)) 
is  the  current  measured  at  the  output  port. 


iMnfflnmiiiifflniiiL 


x  =  0 


jgL 

f"(/)  =  F"V,(r) 


Finite  Element  Mesh 


3  j>-  -4—4  j)  j)  ■ 


Figure  5.  An  external,  transverse  distributed  force  and  generalized 
finite  element  nodal  forces. 


Figure  6.  A  schematic  view  of  the  sensor  port. 


2.3.  Actuator  side 

In  this  work,  the  active  film  under  consideration  is  considered 
to  be  long,  narrow  and  thin,  so  as  to  make  the  length  L\  much 
greater  than  the  maximum  value  of  the  width  b\(x),  which  is 
much  greater  than  the  thickness  hp  (figure  4).  Next,  attention 
is  focused  on  the  case  when  the  active  film  is  driven  by  an 
ac  voltage  V\(t).  The  externally  transverse  distributed  force 
fw(x,  t)  can  be  expressed  in  terms  of  the  drive  voltage  V|(/) 
as 

1  d  2b}(x) 

)  =  -dlxEph'-£^-Vx{  OH(jc),  (15) 

where  H(x)  =  U(x)  —  U(x  —  L\)y  and  U(jc  —  a)  is  the 
Heaviside  step  function.  The  basis  for  (15)  is  detailed  in 
appendix  B. 

In  view  of  ( 1 5),  the  generalized  element  nodal  forces,  ff  (/) 
and  f£(/),  can  be  written  in  terms  of  the  drive  voltage  V\  (/)  as 

ff(0  =®2xl. 

n<')  =  Urf3i£Aj[  ^^N2(jr)d*jv,(»),  0,cq4, 

(16) 

where  N2U)  is  a  shape  function  matrix,  as  discussed  in 
appendix  A.  Further,  QA  =  |  0  ^  x  ^  L\)  denotes  the 

sub-domain  of  [0,  L]  occupied  by  the  sensor  (figure  5).  In  the 
frequency  domain,  ( 1 6)  becomes 

f?(®)=02K|. 

!?(»>  =  ^5^N2C*)drJv,(»).  !!(CQ,, 

(17) 


or 


f»  =  o2xl,  r2(a))  =  r2v]{a>),  si'CtoA.  (i8) 


where  Vi(cu),  f[(cu)  and  (o>)  are  the  Fourier  transforms  of 
Vi(r),  ff(0  and  ^(f),  respectively,  and  Fj  is  the  frequency- 
independent  vector  of  generalized  element  nodal  forces 
defined  as 

„  1  r  d2b,(x) 

FI  =  2^ Jn  2,5!),.  (19) 

By  using  (17)  and  the  assembly  operator,  the  sub- vector  of 
generalized  equivalent  nodal  forces  can  be  written  as  follows: 


rV)  =  a 

e 


^2x1  1 

*2  ) 


V,(w)  =  F>lVI(w), 


eeQA.  (20) 


2.4.  Sensor  side 

On  the  sensor  side,  the  attention  is  focused  on  the  direct 
piezoelectric  effect.  As  described  before,  a  key  to  developing 
a  useful  electromechanical  model  of  a  piezoelectric  model  is 
the  determination  of  changes  in  electrode  charges  when  the 
active  film  is  strained  due  to  the  mechanical  excitation  of  the 
microresonator.  At  the  sensor  port,  it  is  also  considered  that 
the  active  film  is  long,  narrow  and  thin,  so  as  to  make  the 
length  L2  much  greater  than  the  maximum  value  of  the  width 
b2(x ),  which  in  turn  is  much  greater  than  the  thickness  hp 
(figure  6). 

Considering  the  relation  between  the  output  charge  and 
the  output  current  given  by  (B.7)  in  the  frequency  domain,  the 
expression  for  the  output  current  becomes 
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Figure  8.  A  schematic  view  of  the  clamped-clamped  piezoelectric 
microresonator  and  the  finite-element  mesh. 


Figure  7.  Transverse  displacements  of  the  centroidal  axis  and 
generalized  finite-element  nodal  displacements. 


I2(a))  =  j  coQ2(co) 

rL 

=  j  (i)dixEp 


LM 


duo(x,to)  he  d2wo(x 


dx 


ju0(x,  o>)l 

dx2  J 


dr. 

(21) 


The  integral  in  (21)  can  be  split  as 

-U. 

d2W0(X,Ct>)  J  f 


f 

JL- 

f 

Jl-l 


au0(x,w) 

b2(x) - - - dx 

L-Lj  <*X 


9uq(x,  co) 
b2(x) - - - dx, 


b2(x) 

l 

Qe  C  Qs 


ax2 


■dx  ■ 


b2(x) 


dx 

d2Wo(X,  co) 
dx2 


djc, 

(22) 


where  £2$  =  {jc  |  L  —  L2  ^  x  ^  L)  denotes  the  sub-domain 
of  [0,  L]  occupied  by  the  sensor  port  (figure  7).  After  using 
(A20tf)  and  (A20c),  one  obtains 


L. 

L 


3u0(x,  co) 

*,w— 5T-  “ 


-[/, 

■  [/». 


32N2U) 


b2(x) 


dx 2 


j  d'(&>), 
]d! 


(CO). 

(23) 


Hence,  the  expression  for  the  output  current  becomes 
/2(c)  =  [A'(&>)d'(a>)  +  B'(a>)d;(*>)] 

e 

=  Cs(a))d5(aj), 

where  the  sub-matrices  A'(tu)  and  B'(cw)  are  given  by 


(24) 


A'(a>)  =  jWj,  £, 


L 


,  r  x3Ni<*>h 
b2(x) — — — dx. 


dx 


VQ,  C  Qs, 


B'(a>)  =  ~}cojdltEp  J  b2(x)d  ^2(--  dx,  Vfi,  C  Qs. 

(25) 

At  this  point,  it  is  convenient  to  rewrite  the  system  (2)  in  the 
following  block-partitioned  form: 

fA(<u) 


d*(«) 

H  AM(co) 

H'is(a>) ' 

dM(a>) 

»  = 

HM*V) 

HMS(to) 

ds(to)  . 

Hs*(w) 

Hsw(a>) 

Hss((o) 

rM(c v) 
f5(o>) 

(26) 


By  considering  that  in  the  present  case,  both  fM  {co)  and 
fs{co)  are  equal  to  zero,  and  using  (20),  it  is  possible  to  express 
ds(cu)  in  terms  of  the  mechanical  force-displacement  transfer 
function  and  the  driven  voltage  as  follows: 

ds(<w)  =  HM(to)r*(o>)  =  H^MF^co).  (27) 

Finally,  combining  (24)  and  (27),  the  admittance  function 
K21  (co)  relating  input  voltage  to  output  current  can  be  expressed 
as 

Y2](w)  =  =  CV)Hm(»)F\  (28) 

V\(to) 

3.  Experimental  results,  comparisons  and  discussion 

The  developed  coupled  model  enables  one  to  investigate  the 
following:  (i)  the  elastic  stability  of  the  resonator,  (ii)  the 
influence  of  a  constant  axial  force  on  the  transverse  vibrations 
of  a  clamped-clamped  resonator  structure  and  (iii)  the 
influence  of  a  constant  axial  force  on  the  transverse  vibrations 
of  a  free-free  structure. 

3.1.  Pre-stressed  microresonators 

A  case  of  interest  is  one  where  the  geometric  stiffness  is  driven 
by  a  parameter  A;  for  example,  in  the  case  of  a  microresonator 
subjected  to  an  initial  axial  force  Pq  due  to  residual  stresses 
introduced  during  the  fabrication  of  the  resonators  [14-16], 
one  can  write 

Pe  =  \P0,  e=\,2,...,nel  (29) 

where  Pe  is  the  axial  force  in  the  eth  element  and  the  number  of 
elements  is  ne\.  The  geometric  stiffness  is  itself  proportional 
to  this  parameter 

KC=AK^  (30) 

The  eigenvalues  of  the  resulting  problem 

(K* +«$)*,,  =a£M,p.  (31) 

then  become  functions  of  A.  When  com  =  0,  one  gets 
an  eigenproblem  with  eigenvalues  Am  corresponding  to  the 
critical  loads, 

K£Vn  =  (32) 

The  eigenvalue  Aj  of  the  lowest  mode  yields  the  pre-stress 

state 

Pe  er  =  A,  Po.  e=l,2 . rte,  (33) 


where  A  stands  for  ‘actuator’,  M  for  ‘midspan’  and  S  for 
‘sensor’  (figure  8). 


in  which  the  system  buckles.  For  free-free  microresonators, 
one  has  an  unconstrained  structure,  since  the  stiffness  matrix 
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Figure  9.  Variations  of  natural  frequencies  with  an  axial  load  for  an 
AlGaAs  microresonator  (length  =  120  /xm,  width  =  15  /xm, 

L\  =  L2  =  30  /xm,  Q  =  565,  ne]  =  [10  20  10), 

/0  j  =  [1  0427.2.8743,5.6348,9.3146]  kHz). 


Table  1.  Characteristics  of  an  AlGaAs  resonator. 


Layer 

Thickness 
K  0*m) 

Residual  stress 
<70  (MPa) 

Top  AlGaAsiSi 

0.5 

-80 

Middle  AlGaAs  (undoped) 

1.5 

-80 

Bottom  AlGaAs.Si 

2.0 

-80 

is  singular.  In  this  case,  there  are  two  zero  frequencies  for  the 
rigid-body  degrees  of  freedom  and  one  more  zero  frequency 
when  the  compressive  load  equals  the  buckling  load.  To 
remove  the  rigid-body  modes  from  the  solution  of  the  semi- 
definite  eigenvalue  problem  (K£  +  XK *G)<pn  =  one 

can  perform  a  shift  /x  >  0  on  K £  +  XK^  by  calculating 

K£+XK^  +  /xM,  (34) 

and  then,  considering  the  positive  definite  eigenproblem 

(K£  +  XK*g  +  (35) 

where  („  =  <fi„  and  a\  =  +  n  =>  to*  =  a*  -  M- 

3.2.  Influence  of  a  constant  axial  force  on  the  transverse 
vibrations  of  a  clamped-clamped  resonator  structure 

In  figure  9,  for  a  clamped-clamped  AlGaAs  resonator,  the 
variations  of  the  first  four  natural  frequencies  are  shown  with 
respect  to  the  axial  load.  As  expected,  as  the  axial  stretching 
load  increases,  the  natural  frequencies  increase.  Similarly, 
they  decrease  with  the  increase  of  the  compressive  axial 
load. 

When  the  compressive  axial  force  reaches  the  value  of  the 
Euler’s  buckling  load,  PCT  =  1.2489  x  10-2  N,  the  first  natural 
frequency  goes  to  zero.  The  composite  material  properties  and 
layer  thickness  values  for  the  considered  AlGaAs  resonators 
are  listed  in  table  1.  Each  AlGaAs  microresonator  has  a 
particular  orientation  on  the  wafer.  The  orientation  on  the 
wafer  is  indicated,  with  respect  to  a  reference  orientation,  by 
the  so-called  wafer-angle  and  this  in  turn  leads  to  a  certain 
crystallographic  direction  of  the  AlGaAs  structure.  The 
variation  of  the  first  natural  frequency  shown  in  table  2  has 
been  obtained  from  experimental  measurements  for  different 
resonators. 


Figure  10.  Variation  of  the  first  natural  frequency  with  an  axial  load 
P/Pt  for  a  clamped-clamped  AlGaAs  resonator  (Lj  =  L2  = 

25  /xm). 


P/P* 


Figure  11.  An  expanded  plot  of  figure  5  in  the  vicinity  of  the  first 
natural  frequency  of  the  AlGaAs  resonator  predicted  without 
residual  stresses. 


Table  2.  Experimentally  obtained  first  natural  frequencies  for 
different  resonator  lengths  and  wafer  angles. 


Angle 

Length 

15 

21 

27 

33' 

39 

80  /xm 
100  /xm 
120  /xm 

2130  kHz 
1390  kHz 
960  kHz 

2070  kHz 
1360  kHz 
920  kHz 

2040  kHz 
1310  kHz 
900  kHz 

1970  kHz 
1250  kHz 
860  kHz 

1890  kHz 
1200  kHz 
830  kHz 

As  a  representative  example,  a  100  /xm  long  resonator 
with  a  15  /xm  width  is  considered.  From  the  results  shown  in 
figure  10,  it  can  be  seen  that  the  numerically  calculated  value 
of  the  first  natural  frequency  is  1501.5  kHz,  in  the  absence 
of  axial  stresses.  The  first  natural  frequency  value  shifts  to 
1319.3  kHz,  when  the  experimentally  obtained  values  of 
residual  stresses  are  included  in  the  model.  In  figure  11, 
an  expanded  plot  of  a  portion  of  figure  10  is  shown.  The 
horizontal  lines  represent  the  experimentally  obtained  values 
of  the  first  natural  frequency  corresponding  to  the  100  /xm 
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Figure  12.  SEM  of  undercut  introduced  in  a  resonator  anchor  by  the 
etching  processes  during  fabrication  (courtesy,  ARL,  Adelphi,  MD). 


Figure  14.  An  expanded  plot  of  figure  13  around  1400  kHz. 


0 

p/p* 


Figure  13.  Variation  of  the  first  natural  frequency  with  an  axial  load 
P /Pe  for  a  clamped-clamped  AlGaAs  resonator. 

resonator  for  different  wafer  angles.  It  can  be  seen  that 
the  numerically  obtained  values  fall  within  the  range  of 
experimental  measurements,  which  are  listed  in  table  2. 

3.3.  Imperfect  boundary  conditions  due  to  mask 
imperfections  and  fabrication  procedure 

As  a  second  example,  the  effect  of  a  2.5  pun  to  5  pun  undercut 
introduced  in  each  resonator  anchor  by  the  etching  processes 
during  fabrication  [  17]  is  considered  (for  illustrative  purposes 
only,  see  figure  12). 

Resonators  of  105  and  1 10  pun  lengths  were  studied 
to  model  the  effect  of  the  undercut.  In  figure  13,  the 
variation  of  the  first  natural  frequency  with  respect  to  the 
axial  load  is  shown  for  100  pirn,  105  pun  and  110  pun 
long  microresonators.  The  experimentally  obtained  values 
of  the  first  natural  frequencies  are  also  shown  in  the  figure  as 
horizontal  lines  for  making  the  comparisons.  In  figure  14,  an 
expanded  portion  of  figure  13  is  shown  around  the  1400  kHz 
range. 


Here,  the  undercut  has  been  simply  modeled  as  an  increase 
in  the  extension  of  the  microresonator  length.  The  change  in 
the  cross-section  area  is  not  considered  in  this  one-dimensional 
model.  This  effect  has  been  included  in  a  series  of  two- 
dimensional  models  developed  by  the  authors,  and  in  these 
models,  it  is  recognized  that  the  undercut  region  cannot  just 
be  modeled  as  a  beam.  Rather,  a  more  refined  finite-element 
model  with  plate  elements  needs  to  be  used  to  model  more 
realistic  boundary  conditions  which  include  changes  in  the 
cross-section  area. 

3.4.  Material  property  variations  from  manufacturing 
uncertainties 

As  a  third  example,  the  effect  of  a  ±10%  varying  elastic 
modulus  across  the  length  of  the  multi-layered  microresonator 
is  considered.  For  micromechanical  resonators,  these 
uncertainties  in  material  properties  come,  in  large,  partly  due 
to  the  manufacturing  process  [17].  A  resonator  of  120  pun 
length  is  studied  to  model  the  effect  of  the  varying  elastic 
modulus.  In  figure  15,  the  variation  of  the  first  natural 
frequency  with  respect  to  the  axial  load  is  shown  for  90 %E,  E 
and  1 10%£.  The  experimentally  obtained  values  of  the  first 
natural  frequencies  are  also  shown  in  the  figure  as  horizontal 
lines  for  making  the  comparisons.  In  figure  16,  an  expanded 
portion  of  figure  15  is  shown  around  the  900  kHz  range.  For 
all  wafer  angles,  the  measured  average  residual  stresses  are 
the  same  and  they  equal  -80  MPa. 

The  residual  stresses  are  based  on  lattice  structure 
measurements,  and  these  values  were  provided  to  the  authors 
by  researchers  in  the  Maryland  MEMS  Laboratory.  Wafer 
bow  measurements  have  also  been  used  to  measure  residual 
stresses  in  this  type  of  multi-layered  microstructures  [18]. 

3  J.  Predictions  of  frequency  response 

The  shapes  of  the  drive  and  sense  electrodes  affect  the 
resonator  admittance  through  £|(jc)  and  bi(x)  in  (19)  and 
(25).  For  the  resonators  shown  in  this  work,  maximum 
electromechanical  coupling  is  desired;  that  is,  the  electrodes 
must  be  shaped  such  that  Y2\(co)  is  maximized.  This  may  be 
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Figure  15.  Variation  of  the  first  natural  frequency  with  an  axial  load 
P/Pc  for  a  clamped-clamped  AlGaAs  resonator. 


Figure  16.  An  expanded  plot  of  figure  15  around  900  kHz. 


Figure  17.  The  response  of  a  400  /xm  AlGaAs  resonator  for  three 
different  values  of  the  axial  load  (Q  =  565,  Pc  =  7.4930  x 
KT4  N.  L,  =L 2  =  100  Mm,  NE  =  [10  20  10],/,  = 

[81.59  93.84  104.55]  kHz). 


^  -  20  fim 


Figure  18.  PZT  microresonator  characteristics. 


Table  3.  Material  properties  for  finite-element  modeling. 


Material 

Young’s 

modulus 

E  (GPa) 

Density 
p(kg  m-3) 

Poisson’s 

ratio 

V 

Pt 

160-190 

18000 

0.35 

PZT 

25  or  35 

8  800 

0.30 

Si02 

100 

2200 

0.27 

achieved  by  clipping  the  electrodes  at  the  quarter  beam  points, 
as  depicted  in  figure  1(6).  For  this  electrode  geometry,  b\{x) 
and  62(x)  can  be  written  as 


b\(x)  =  b  |^U(x)  —  U 


-  U(jr  -  Z,)l  =  bg2(x). 


(36m) 

(366) 


where  b  is  the  nominal  width  of  both  electrodes  and  U(x)  is 
the  Heaviside  unit  step  function.  By  using  (36),  the  terms  in 
(19)  and  (25)  involving  the  electrode  shapes  simplify  to 

F5  =  Qrfji  Eph'b )  jf  ^£^N2(jt)  dx  (37) 
and 


/ 1  \  2  f  9N,(x) 

\e(co)  =  jo;  y-hed3lEpbJ  —  J  g2(x)  ^  dx, 

B'(a))  = -jco(^heditEpbj  J  gz(x)- -  dr. 


(38) 


In  figure  17,  for  a  400  /xm  doubly  clamped  AlGaAs 
resonator  with  quarter-beam  electrodes,  the  response 
(magnitude  and  phase  of  the  normalized  admittance  function 
Y2\{(*>) / (\di\ Eph'b)1)  is  shown  for  three  different  values 
of  the  axial  load.  As  expected,  as  the  axial  stretching  load 
increases,  the  first  natural  frequency  increases.  Similarly,  it 
decreases  with  the  increase  of  the  compressive  axial  load.  The 
composite  material  properties  and  layer  thickness  values  for 
the  considered  AlGaAs  resonators  are  listed  in  table  1. 

As  a  second  example,  a  200  /xm  clamped-clamped  PZT 
resonator  with  quarter-beam  electrodes  (see  figure  18)  is 
considered  and  the  responses  obtained  for  three  different 
values  of  the  axial  load  are  shown  in  figure  19.  Again, 
as  expected,  as  the  axial  stretching  load  increases,  the  first 
natural  frequency  increases.  Similarly,  it  decreases  with  the 
increase  of  the  compressive  axial  load.  The  material  properties 
and  geometric  properties  used  for  finite-element  modeling  of 
the  considered  PZT  resonator  are  listed  in  tables  3  and  4, 
respectively. 
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Appendix  A.  The  finite-element  model 


Frequency  (kHz) 


Figure  19.  The  response  of  a  200  fim  PZT  microresonator  for  three 
different  values  of  the  axial  load  ( Q  =  4000,  Pe  =  3.3338  x  10"3  N, 
Lx  =  L2  =  50  Mm,  NE  =  [10  20  101,  /i  =  [255  303  344]  kHz). 


Table  4.  Geometric  properties  for  finite-element  modeling. 


Layer 

Thickness  /ik 

Si02 

1.0602  Mm 

Bottom  Pt 

135.0  nm 

PZT 

0.5303  Atm 

Top  Pt 

200.0  nm 

4.  Summary 

In  this  work,  a  semi-analytical  computational  mechanics 
model  of  a  composite  microstructure  with  piezoelectric 
actuation  and  piezoelectric  sensing  has  been  developed 
as  a  design  tool  for  microresonators.  The  developed 
dynamic  model  of  microresonators  accounts  for  structural 
properties  and  the  electromechanical  coupling  effect  through 
finite-element  analysis.  The  dynamic  admittance  model  is 
derived  by  combining  the  linear  piezoelectric  constitutive 
equations  with  the  modal  transfer  function  of  the  multi-layered 
microresonator  structure.  The  resonator  receptance  matrix  is 
constructed  through  modal  summation  by  considering  only  a 
limited  number  of  dominant  modes.  The  electromechanical 
coupling  determination  at  the  input  and  output  ports  makes  use 
of  the  converse  and  direct  piezoelectric  effects.  The  developed 
model  has  been  validated  by  comparing  it  with  results 
available  in  the  literature  for  clamped-clamped  resonators. 
The  numerical  results  are  found  to  be  in  good  agreement  with 
the  experimental  measurements.  The  numerical  simulations 
show  that  the  consideration  of  axial  loads  is  important  to 
predict  the  natural  frequencies  of  the  resonators  studied 
in  the  experiments.  A  detailed  discussion  of  modeling 
considerations  has  also  been  presented.  The  microresonators 
studied  in  this  work,  which  are  used  as  micromechanical  filters, 
are  important  for  mobile  communication  systems  and  signal 
processing  applications.  It  is  believed  that  the  numerical  and 
analytical  efforts  presented  in  this  work  can  be  used  as  a  basis 
to  develop  design  tools  for  such  systems. 


In  this  appendix,  a  geometrically  nonlinear  finite-element 
analysis  (e.g.  [19,  20])  is  carried  out  to  study  transverse 
vibrations  of  clamped-clamped  and  free-free  resonators 
subjected  to  constant  axial  loads.  This  is  done  in  the  context 
of  a  planar  beam,  but  a  similar  development  can  be  carried 
out  for  the  three-dimensional  case.  As  discussed  earlier, 
the  consideration  of  axial  loads  is  important  to  predict  the 
first  natural  frequencies  of  the  resonators  observed  in  the 
experiments.  In  addition,  as  also  shown,  the  developed  model 
and  the  numerical  implementation  can  be  used  to  understand 
the  influence  of  uncertainties  associated  with  the  fabrication 
process.  The  effect  of  curvature  shortening  (e.g.  [21])  is 
incorporated  in  the  developed  numerical  code. 


A.l .  Kinematics  of  the  motion 


The  equations  governing  the  bending  of  Euler-Bemoulli 
beams  with  moderately  large  rotations  but  with  small  strains 
can  be  derived  using  the  following  linearized  displacement 
field 

dw0(x,  t) 

ui=u0(xj)-z - - - ,  U  2  =  0,  u3  =  w0(x,t), 

dx 

(A.l) 


where  (u\,U2,u3)  are  the  total  displacements  along  the 
coordinate  directions  (jc,  y,  z),  and  u0  and  u>o  denote  the  axial 
and  transverse  displacements  of  the  centroidal  (neutral)  axis 
(y  =  z  =  0).  Equation  (A.l)  is  used  in  the  majority  of 
existing  nonlinear  formulations.  Although,  linearization  of 
beam  kinematics  is  apparently  incongruent  with  the  notion 
of  geometric  nonlinearity,  the  linearized  field  provides  an 
adequate  description  for  incremental  static  and  dynamic 
solutions.  In  a  geometrically  nonlinear  analysis,  using  the 
finite-element  method  formulation,  it  is  customary  to  use  a 
nonlinear  strain-displacement  relation  of  the  form  (sum  on 
repeated  subscripts  is  implied) 


Due  to  the  assumption  of  small  strains,  no  distinction  is  made 
here  between  the  material  and  the  spatial  coordinates  [22]. 
In  the  case  of  planar  motions,  substitution  of  the  linearized 
displacement  field  (A.l)  into  (A. 2)  yields 


£n 


—  2z 


duo  d2Wo 
dx  dx2 


z 


d2w0 

HI*"* 


(A. 3) 


and  all  other  strains  are  zero.  Note  that  the  notation  x\  =  jc, 
*2  =  y  and  JC3  =  z  is  used.  The  first  term  in  (A. 3)  represents 
the  centroidal  strain  of  the  element  according  to  the  Lagrangian 
(Cauchy-Green)  strain  tensor,  and  the  last  term  represents  the 
additional  axial  strain  due  to  flexure  for  fibers  at  a  distance 
y  from  the  centroid  of  the  cross-section.  For  the  considered 
microresonators,  the  longitudinal  inertia  is  small  compared 
with  the  restoring  force,  and  hence,  based  on  this  hypothesis 
it  follows  that  uq  =  0(wl)  (i.e.  u0  is  of  the  order  of  Wq)  [23]. 
Making  use  of  this  assumption  and  neglecting  terms  of  order 
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Figure  20.  Geometry  of  a  typical  finite  element. 


higher  than  2;  that  is,  ignoring  terms  that  involve  the  square  of 
dup/dx  and  retaining  the  square  of  dwp/dx  (which  represents 
the  rotation  of  a  transverse  normal  line  in  the  beam)  yields 


These  strains  are  known  as  von  Karm&n  strains. 


(A. 4) 

(A. 5) 


A. 2.  Strain  energy 


For  a  given  element,  the  strain  energy  Ue  is  given  by 

U‘(t)=  [  [  Ojj  dEjj  dV . 

JV'  Jc^ 


(A. 6) 


Here  Ve  is  the  volume  of  the  element,  <ty(jr,z,  t)  and 
Cij(x,  z,t)  denote  the  Cartesian  components  of  the  stress  and 
the  Green  strain  tensors,  respectively.  Due  to  the  assumption 
of  small  strains,  no  distinction  is  made  here  between  the 
Cauchy  and  the  second  Piola-Kirchhoff  stress  tensors  [22]. 
The  total  volume  of  the  laminated  element  is  given  by  the 
tensor  product  Ae  x  Qet  where  Ae  is  the  cross-sectional  area, 
and  Cle  =  [xe,  jc,+i]  defines  the  length  of  the  element.  Here, 
xe  and  xe+\  are  the  jc-coordinates  of  the  left  and  right  nodes  of 
the  element,  as  shown  in  figure  20. 

In  view  of  the  explicit  nature  of  the  assumed  displacement 
field  (A.l)  in  the  thickness  coordinate  z  and  its  independence 
of  coordinate  y,  the  volume  integral  can  be  expressed  as  the 
(tensor)  product  of  integrals  over  the  length,  xe  ^  x  ^ 
and  the  cross-sectional  area  of  the  element: 


[  (*)  d  V 
Jve 


-CL 


(  )d A  Ax. 


Therefore,  the  expression  for  the  strain  energy  can  be 
simplified  as  follows  (for  a  beam  element,  the  only  nonzero 
components  of  strain  and  stress  are  e n  =  exx  and  a\\  =  (Txx): 


u'(t)= L  L  l  (ff"d£«)<Mdt’  (a-7) 

where  axx  is  the  axial  stress  and  exx  is  the  axial  strain  given  by 
(A. 4).  Assuming  that  each  layer  is  isotropic  with  respect  to  its 
material  symmetry  lines  and  linearly  elastic,  the  axial  stress 
and  strain  are  related  by  Hooke’s  law 


t )  =  0,  *=1.2 . ft  layers 

(A. 8) 


where  ft|ayere  denotes  the  total  number  of  layers  and  £<*>( x)  is 
the  modulus  of  elasticity  of  the  ki h  layer.  After  substituting 
(A. 8)  into  (A. 7)  and  performing  integration  over  the  strain 
field  f.xx>  the  result  is 


•»-/  /  / 

Jnr  J  a,  Jc, 

-\U, 


(£e„  de„)  dA  dr 
Eel,  d A  dx. 


(A.9) 


After  substituting  (A.4)  into  (A.9),  the  expression  obtained  is 


U&P 


du0  d2w0 

\ - —z 

dx  dx2 


d2Wp  /  dw0\2  t  du  0  /  dll)0\2 

dx2  \  dx  )  Z  dx  \  dx  J 


dA  dx. 


(A.  10) 


The  higher-order  term  1/4 (dw0/dx)4  can  be  neglected  in  the 
above  expression.  Although  the  strains  are  continuous  through 
the  thickness,  stresses  are  not,  due  to  the  change  in  material 
coefficients  through  the  thickness  (i.e.  each  lamina).  Hence, 
the  integration  of  stresses  through  the  laminate  thickness 
requires  laminawise  integration.  Integrating  (laminawise)  over 
the  cross-sectional  area  Ae  and  noting  that  since  z  is  measured 
from  the  neutral  axis  of  the  multi-layered  cross-section,  all 
integrals  of  the  form  f  zdA  must  vanish,  one  obtains 


where  (EA)eff(jc)  and  ( EI)en(x )  denote  the  effective 
extensional  and  bending  stiffness  of  a  typical  (generic)  finite 
element,  respectively.  One  can  note  that  the  first  two  integrals 
in  (A.  11)  represent  the  linear  strain  energy  while  the  third 
integral  is  the  contribution  from  the  nonlinear  component  of  the 
strain,  l/2(dw0/dx)2.  In  order  to  derive  the  stiffness  matrix 
from  the  expression  for  the  strain  energy  given  by  (A.  11), 
it  is  necessary  to  express  the  displacement  fields  uo(x ,  /)  and 
wo(x,  t )  in  terms  of  the  generalized  nodal  displacements  d‘(t), 

i  =  1,2 . 6  (see  figure  21).  This  can  be  accomplished  by 

assuming  a  displacement  field  for  up(x ,  t )  and  wo(*.  /).  By 
letting 


uo(x ,  t)  =  ap{t)  +ai(/)*,  (A.  12) 

W0(x,t)  =  bo(t)  +  b,(t)x  +  bi{t)x2  +  b3(t)x3,  (A.  13) 

and  using  the  boundary  conditions, 
uo(x.  01,=,,  =  d\( t), 
uto(x, /)!,=,,  =d'2(t), 

wq(x,  /)!,=,„,  =  d'5(t), 


uo(x,  01,=,,.,=  d^(t),  (A. 14) 

a 

—  iuo(x.  01,.,,  =  </3(0, 

OX 

—  U)0(x, /)!,=,„,  =<(/). 

(A.  15) 
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Figure  21.  Generalized  nodal  displacements  and  finite-element 
basis  functions. 


it  can  be  shown  that 

u0(x,t)  =  [Nl(x)  N<(x)] 

where 

Ni(x)=  1-f, 


<(/) 

did) 


N*(x)  =  - 

'e 


w0(x,t)  =  lN2(x)  N}(x)  N5(x)  N6(x)) 


d'2«) 

did) 

d'5d) 

did) 


where 


N2(x)=  1  -3-J-+2-5-, 

le  *g 

X2  X3 

N,(x)  =  37T-27T, 


N3(x)  =  x r-2-  +  -, 

W,)  =  -77  +  7r 


and 


uo(x,t)  =  Ni(ar)d^f) 

w0(x,t)  =  N2(jr)d2(/), 


dx 


dx 


dx2 


3N'('k(0. 

dx 

(A.20  a) 

(A. 20b) 

(A.20  c) 

In  view  of  (A.  18)-(A.20),  the  element  strain  energy  Ue  can  be 
written  as 


u'm=iL‘' 

♦u 


3N[U)aN|(^)  . 

" lTd)EA '  »l  .u..'df(,)d» 
dx  dx 


A‘‘  (t)Elr 


3N[WaN2(ar) 


dx 


dx 


dUOdx 


[- 


"KwTTw]'-. 


(A. 21) 


^N[(jc)  3N2(jr)Jtf/ 
ax  ax 

It  is  to  be  noted  that  even  for  relatively  large  deflections,  the 
term  EA,/lrdN](x)/dxd'(t)  may  be  treated  as  a  constant 
equal  to  the  axial  tensile  force  in  the  element.  Hence, 
introducing 


N;w  r  i  ilk 

h  dx  ,V  'L  ‘e  l A  }< 


did) 

did) 


(A.  16) 


=  ^[dSd)-di(,)], 

le 

it  is  obtained  that 


(A. 22) 


r(0 


are  the  linear  Lagrange  basis  functions,  le  =  xe+\  —  xe  is  the 
length  of  the  element  ft,  and 


[/,"'• 

•U 


aNr  (x)  aN2(x) 


dx 


dx 


+  -d  ‘/(OP. 


A. 3.  Kinetic  energy 


aN[w  aN2u) 


dx 


dx 


djrjd^(r) 

■*  j  d2(0- 


(A.23) 


In  this  work,  the  effect  of  rotary  inertia  is  ignored.  Consistent 
with  this  kinematic  assumption,  the  kinetic  energy  is  given  by 


(A. 17) 


~dw0(x,  t ) 
~dt 


dx, 
(A. 24) 

where  meff  (x)  is  the  effective  mass  per  unit  length  of  the  multi¬ 
layered  element.  In  view  of 


are  the  cubic  Hermite  basis  functions.  Symbolically,  (A.  16) 
and  (A.  17)  can  be  written  respectively  as 


^,o=NiW®o=NiWd;(a 


dt 

dw0(x,t) 


=  N2(x) 


dt 

dd‘2(t) 


=  N2(x)d‘(f), 


(A. 18) 


(A. 19) 


where  Nj(x)  and  N2(x)  are  referred  to  as  shape  function 
matrices.  Furthermore, 


dt  '  dt 

the  element  kinetic  energy  Te  can  be  written  as 

T  .  , 

1 


(A. 25a) 
(A.25b) 


r‘d)  =  X-A\t (f)  jjf  m^(jc)N[(ac)N,(jt)<ujd^(0 


[I, 


meirC*)N2  (jt)N2(jt)  dx 


]d'2(0. 


(A. 26) 


A. 4.  Virtual  work  due  to  nonconservative  distributed  forces 

The  virtual  work  done  by  the  nonconservative  distributed 
forces  has  the  form 

<5^nr(/)=  [  fu(x,t)8u0(x,  t)  dx  +  [  fw(x,t)8w0(xyt)dx, 
Jnr  Juf 

(A. 27) 
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where  fu(x ,  t )  is  the  externally  applied  axial  force  per  unit 
length,  Suq(x ,  t)  is  an  arbitrary  axial  virtual  displacement, 
fw(x,t)  is  the  externally  applied  transverse  force  per 
unit  length  and  fiw0(x,  t )  is  an  arbitrary  transverse  virtual 
displacement.  In  view  of  (A.  18)  and  (A.  19),  the  virtual 
displacements  &uq(x,  t)  and  &w0(x,  t )  can  be  written  as 


&uo(x,  t)  =  N,(jc)ad*(/) 

(A. 28) 

and 

&w0(x,  t )  =  N2(jc)<5d2(f), 

(A. 29) 

respectively.  On  using  (A. 28)  and  (A. 29)  in  (A. 27),  the  virtual 
work  due  to  the  nonconservative  distributed  forces  has  the  form 

sw"c(i)=\L 

+  |^  fw(x, ON2C*) dx | £dj(f) 

(A. 30) 

where  ff  (f)  and  fj(f)  are  the  generalized  element  nodal  forces. 
Note  that  the  integration  domain  above  is  the  small  sub-domain 
(0  ^  x  ^  le)  of  a  typical  finite  element. 


the  geometrical  stiffness  K^.  The  elastic  stiffness  matrix 
is  the  same  as  that  used  in  the  linear  analysis  and  its  nonzero 
sub-matrices  are  given  by 


Ke„  =  [  (^)effW 
Jn, 


dNi(jc)  9N[ (jc) 


dx 


dx 


dx 


and 


K22  =  [ 
Jn, 


,  x3N[(x)3N2U)  j 
(E/)ef {(x) — r - - - dx. 


dx 


dx 


(A. 37) 


(A. 38) 


The  nonzero  component,  KeG22>  °f  the  geometrical  stiffness 
matrix,  K^,  is  given  by 


3N[U)  3N2(-t) 
dx  dx 


dx. 


(A. 39) 


It  is  important  to  mention  that  the  geometrical  stiffness  matrix 
either  increases  or  decreases  the  direct  stiffness  coefficients 
(diagonal  terms),  depending  on  the  sign  of  the  axial  force  Pe. 
After  substituting  (A. 32),  (A. 33)  and  (A. 36)  into  (A. 31),  one 
obtains 


A. 5.  Equations  of  motion  in  the  time  domain 


M'd'(r)  +  [K'£  +  Kg]d'(/)  =  f(0,  (A. 40) 


In  this  section,  Lagrange’s  equations  [  1 3,  24]  are  used  to  derive 
the  equations  of  motion.  Treating  the  nodal  displacements 
d'r(o  =  {df  (O.d^CO}  as  ^e  generalized  coordinates, 
Lagrange’s  equations  take  the  form 


d  /  ar^\  <W_ 

dt  VacK  /  ad77  +  dd'T 

where 


(A. 31) 


rr(o  =  *2r(')|- 


After  using  (A. 26),  one  can  write 

_d  /  37^\  _  I’M',  02x4l  |  dj(0 1 
dr  \9#t)  ~  [o4x2  MJ2J  J 


M'd'(r), 


(A. 32) 


(A. 33) 


where  is  the  element  consistent  mass  matrix,  and  its 
nonzero  sub-matrices  and  M22  are  defined  as 


and 


M'u=  [  OTc(r(jr)N[(Ar)N1(Ar)dx 

Jn, 

M'22=  f  me(f(x)N[(jc)N2(jr)<Lt. 

Jn, 


(A.34) 


(A. 35) 


In  this  work,  since  the  system  is  a  natural  system  (the  kinetic 
energy  has  the  quadratic  form  shown  in  (A. 26))  [20,  21] 
and  its  mass  distribution  does  not  depend  on  the  generalized 
coordinates;  hence  dTe/ddeT  =  06xj.  By  using  (A. 23),  the 
last  term  on  the  left-hand  side  of  (31)  can  be  written  in  the 
following  form: 


dUe 

ddeT 


\\*sn  ®2*4' 

j  04x2  ^E22 

[K‘£+K'c]d'(0. 


02x2 

04x2 


02x4  "j  |  (d '(f)' 

K'C22  Jtid5(0 

(A. 36) 


Thus,  one  can  clearly  see  that  the  total  stiffness  of  the 
element  consists  of  two  parts,  the  elastic  stiffness  and 


which  is  the  equation  of  motion  for  element  e  in  the  global 
inertial  reference  system.  Expanding  this  equation  of  motion 
to  system  size  and  combining  all  the  element  equations,  the 
equation  of  motion  for  the  entire  system  can  be  written  in  the 
following  form: 

Md(r)  +  Cd(r)  +  Kd(r)  =  f(r),  (A.41) 


where 


K  =  K£  +  Kc  (A. 42) 

is  the  global  stiffness  matrix,  and  M,  Kf,  and  f(r)  are 
the  (global)  mass  matrix,  elastic  stiffness  matrix,  geometrical 
stiffness  matrix  and  force  vector,  respectively.  All  these 
matrices  and  the  force  vector  of  the  complete  system  are 
obtained  from  the  corresponding  individual  element  matrices 
and  load  vector  by  using  the  direct  stiffness  method;  that  is, 

M  =  A  (MO.  K £  =  A  (K£),  Kc  =  A  (K'c) 

e=l  e=  1  r=l 


and 


f(0  =  a  (raw. 

e=\ 

where  A  is  the  so-called  assembly  operator  [25]  and  nc\  is  the 
number  of  elements  in  our  model. 

Here,  the  Rayleigh  proportional  damping  is  introduced  in 
the  equation  of  motion  with  the  following  formulation  for  the 
damping  matrix: 


C  =  aM  +  /3K.  (A. 43) 

In  (A. 43),  a  and  fi  are  coefficients  which  can  be  obtained  from 
the  modal  properties  of  the  microresonator. 
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Appendix  B.  Actuator  and  sensor  models 

The  derivation  of  the  electromechanical  coupling  at  the  input 
and  output  ports  makes  use  of  the  constitutive  equations 
of  a  linear  piezoelectric  material.  The  length  direction  of 
the  piezoelectric  film  is  chosen  along  the  x-axis,  the  width 
direction  is  chosen  along  the  y-axis  and  the  height  direction  is 
chosen  along  the  z-axis.  This  choice  of  variables  implies 
that  the  structure  of  the  constitutive  equations  of  a  linear 
piezoelectric  material  in  the  ‘T-E’  form  is  [26] 


+  J31E3,  (B.la) 

D3  =  ^31  T\  +  £33  E$.  (BAb) 


Here,  the  contracted  (Voigt)  notation  is  used  to  simplify  the 
presentation  of  second-  and  fourth-order  tensors.  In  (B 1 )  s 
Exx  and  T\  =  axx  are  the  strain  and  stress  in  the  jr-direction 
of  the  beam;  O3  and  £3  are  the  electrical  displacement  and 
electrical  field  along  the  piezoelectric  film  in  the  z-direction 
and  sfj,  d^\  and  are  the  x-axis  mechanical  compliance  at 
a  constant  electric  field,  the  transverse-piezoelectric-coupling 
coefficient  and  the  permittivity  at  a  constant  stress  of  the  film, 
respectively. 

At  the  input  port  (i.e.  the  actuator  side),  through  the 
converse  piezoelectric  effect,  the  voltage  applied  across  the 
active  film  causes  a  forced  strain  equal  to 

5,(0  =  *i£3(0.  (B.2) 


where  £3(0  is  the  electric  field  across  the  active  film,  which 
is  given  by 


£3(0  = 


v,(0 


(B.3) 


This  forced  strain  is  caused  by  the  electric  field,  which  is 
uniform  through  the  active  film;  hence,  the  forced  strain  is 
also  uniform  in  the  active  film.  This  means  that,  if  we  restrain 
the  active  layer  from  an  extension,  the  forced  stress  is  also 
uniform  through  the  layer,  and  it  is  equal  to 


r,(o  =  ^  =  st(t)Ep, 


(B.4) 


where  Ep  is  Young’s  modulus  of  the  piezoelectric  film.  Since 
the  film’s  offset  from  the  neutral  axis,  the  forced  stress  given  by 
(B.4)  is  translated  into  a  distributed  internal  bending  moment 
To  obtain  an  approximation  for  the  distributed 
moment,  one  can  follow  along  similar  lines  of  earlier  work. 
The  distributed  bending  moment  is  given  by  the  product  of  the 
lateral  force  caused  by  the  piezoelectric  film,  T\(t)hpb\{x), 
and  the  distance  between  the  film  and  the  neutral  axis  of 
the  microresonator.  Assuming  that  the  thickness  of  the 
piezoelectric  layer  {hp)  is  relatively  small  compared  to  the 
thickness  of  the  resonator  substrate  (he),  after  using  (B.2)- 
(B.4),  the  expression  for  M (jc,  t)  becomes 

A#(*.i)  =  ffiWMitoly  =  [5,(0£,M,<*)]y 


(B.5) 


It  is  known  that  an  externally  transverse  distributed  force, 
fw(x>  0*  leads  to  an  internal  bending  moment  at  each  section 
of  the  microresonator  structure,  but  here  one  has  the  reverse; 
one  has  an  internally  generated  distributed  moment  inside  the 


microresonator  structure,  and  it  is  needed  to  bring  that  outside 
as  an  externally  distributed  force  that  is  structurally  equivalent 
to  M ( x ,  /).  From  [27],  it  follows  that  the  relationship  between 
M (jc,  t)  and  fw( jc,  t )  is  given  by 

d2M(x,t) 

fw(x,t)=  -dx2~-  (B.6) 

At  the  output  port  (i.e.  the  sensor  side),  the  output  current, 
/2(f)*  can  be  written  in  terms  of  the  output  charge,  Q2(t),  as 
follows: 

/2W-XC2O). 


d  / 

The  output  charge  is  found  by  integrating 
displacement,  £>3(jc,y,f).  over  the  area  of 
electrode: 

rh(*) 


the 

the 


(B.7) 

electric 

sensing 


Qi(t) 


-u 

- LJ 


£>3  dy  dx 


(dy\  Ti  +  £33  £3)  dy  dx. 


(B.8) 


In  the  case  of  a  noncollocated  sensor,  the  voltage  across  the 
sensor  is  maintained  at  £3(0  =  0.  After  using  this  fact  and 
appealing  again  to  the  actuator  equation  (B.8)  simplifies  to 

rhix) 


02(0 


-f  f 

J  L—Li  Jo 


d3]T\(x,t)dydx. 


(B.9) 


Here,  the  stress  in  the  sensing  layer,  T\(x,t),  can  be  written 
in  terms  of  the  longitudinal  strain,  5,(x,  /),  and  the  Young’s 
modulus  of  the  piezoelectric  film,  Ep ,  as 

T,(jc,0  =  Ep5\  (jc,  t).  (B.10) 

Since  we  assume  that  the  thickness  of  the  piezoelectric  sensing 
layer  ( hp )  is  relatively  small  compared  to  the  thickness  of  the 
resonator  substrate  ( he ),  the  longitudinal  strain  in  this  layer  is 
approximately  equal  to  the  longitudinal  strain  at  the  surface  of 
the  microresonator  structure.  This  total  strain  consists  of  both 
a  bending  strain  and  an  extensional  strain: 

he  d2w0(x ,  /) 


5®end(jt,0  =  — 


Sf-(x.0«2^  + 


dx 


dx 2  ’ 

1  /  dw0(x, 

2  V  dx~ 


(B.lla) 


(B.ll  b) 


In  this  work,  the  contribution  from  the  nonlinear  component 
of  the  extensional  strain  to  the  output  charge  is  neglected. 
When  the  input  at  the  drive  electrode  has  a  dc  bias,  in  addition 
to  the  time  varying  excitation,  the  effect  of  the  nonlinear  term 
is  included  in  the  model  as  a  constant  axial  load.  With  this 
assumption,  the  output  charge  is  found  by  combining  (B.9)- 
(B.l  1),  resulting  in 

rM*>  ra,,-  hrd2w0' 


02(0  =  d)i  E . 


=  dji  Ei 


rL  r^pMo  ht  32ui01 
Jl-lJo  L  Ojt  2  »jr*J 

OH 


du0(x,  t ) 
dx 


Dx*  'dytbr 
he  d2wp(x,  0 


dx2 


dx, 


and,  in  the  frequency  domain, 
rL 


Qi{oj)  =  d^\Ep  f  b2(x) 
Jl-Li 

«[ 


l-l7 

[duo(x,co)  hed2wo(x,o)) 


]*■ 


(B.12) 


(B.13) 


ax  2  ax1 

where  w0(x,a>)  and  u0(x,w)  are  the  Fourier 

transforms  of  02(f).  wo(x,  l)  and  uq(x,  i),  respectively. 
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Abstract 

Piezoelectrically  actuated  and  sensed  microscale  resonators  have  been 
widely  studied  for  filtering,  communication  and  sensing  applications.  In  this 
effort,  to  characterize  these  resonators,  the  nonlinear  frequency-response 
behavior  of  these  resonators  is  examined  and  parametric  identification  is 
carried  out.  A  nonlinear  beam  model  is  used  with  a  single-mode 
approximation  to  produce  a  forced  Duffing  oscillator.  Nonlinear  analysis  is 
used  to  obtain  the  frequency-response  equation,  and  this  equation  is  used 
along  with  a  least-squares  minimization  scheme  to  identify  the  linear  and 
nonlinear  parameter  values  in  oscillator  models  describing  the  microscale 
structures.  A  linearized  analytical  model  of  the  stepwise  axially  varying 
resonator  is  also  used  to  obtain  additional  system  parameters.  The 
experimentally  identified  parameter  values  are  found  to  be  in  agreement  with 
predicted  values  obtained  from  a  nonlinear  beam  model.  Parameter  values 
obtained  from  multiple  sets  of  data  for  PZT  and  AlGaAs  microresonators 
are  used  to  observe  trends  with  respect  to  a  variety  of  operating  conditions. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Piezoelectric  microscale  resonators  and  resonator  arrays 
are  currently  being  developed  through  a  number  of 
different  efforts.  These  devices  are  manufactured  by 
using  microelectromechanical  systems  (MEMS)  fabrication 
techniques  with  the  goal  of  replacing  the  bulky  macroscale 
resonators  currently  in  use.  With  microscale  dimensions,  they 
require  only  a  fraction  of  the  space,  and  in  addition,  the  power 
requirements  are  low  compared  to  those  of  their  macroscale 
counterparts.  A  number  of  different  resonator  devices  of 
various  geometries  are  currently  being  developed.  A  few 
of  these  include  clamped-clamped  resonators  [1],  film  bulk 
acoustic  wave  resonators  (FBARs)  [2],  ring-shaped  contour¬ 
mode  resonators  [3]  and  various  disk  resonators  [4,  5].  During 
their  development,  the  presence  of  nonlinear  behavior  has  been 
observed  in  many  of  these  devices  [6,  7].  Modeling  and 
identification  of  the  nonlinear  behavior  of  dynamic  systems 
has  been  the  focus  of  a  large  number  of  studies.  These 
efforts  include  both  non-parametric  identification  methods 
[8]  and  parametric  identification  methods  [9-11].  Here,  a 
parametric  identification  scheme  is  developed  and  applied  to 
study  the  nonlinear  dynamic  behavior  of  clamped-clamped 
microresonators. 


Ayela  and  Fournier  [12]  proposed  a  method,  in  which  a 
Duffing  oscillator  model  is  used  to  determine  the  parameter 
values  for  microstructures  exhibiting  nonlinear  behavior. 
In  their  identification  scheme,  seven  parameter  values  are 
determined  by  using  seven  equations.  These  equations 
correspond  to  the  amplitude  and  frequency  of  the  maxima 
associated  with  the  increasing  (or  forward)  and  decreasing 
(or  backward)  frequency  sweeps,  the  frequency  difference 
between  the  location  of  the  maximum  determined  during  the 
increasing  or  forward  frequency  sweep  and  the  point  where  the 
response  drops  to  the  lower  branch  of  the  response  curve  and 
the  critical  force  and  amplitude  values.  These  critical  values 
correspond  to  the  separation  of  the  responses  observed  during 
increasing  and  decreasing  frequency  sweeps  into  two  unique 
curves.  By  using  a  large  number  of  data  points,  the  method 
proposed  in  this  paper  allows  for  a  better  representation  of 
the  response  in  the  presence  of  noise  and/or  small  deviations 
from  the  ‘ideal’  case,  where  the  experimental  data  align 
perfectly  with  a  frequency-response  curve  determined  by  using 
the  forced  Duffing  equation.  Furthermore,  the  scheme  does 
not  require  data  points  collected  during  the  backward  (or 
decreasing)  sweep  of  the  excitation  frequency.  Also,  the 
method  developed  here  provides  additional  information  such 
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Figure  1.  Structure  of  the  piezoelectric,  microscale 
clamped-clamped  resonator. 


h-  I- 


Figure  2.  PZT  resonator  and  geometry. 

as  the  average  residual  stress  levels  in  a  microscale  structure. 
A  parametric  identification  scheme  was  also  proposed  by 
Malatkar  and  Nayfeh  [  1 3]  to  study  systems  exhibiting  similar 
nonlinear  dynamic  behavior.  However,  unlike  the  scheme 
described  in  this  paper,  the  identification  scheme  of  Malatkar 
and  Nayfeh  was  meant  to  be  applied  to  frequency-response 
data  that  did  not  exhibit  jumps. 

In  the  next  section,  the  microscale  structures  of  interest 
are  discussed  and  the  experimental  arrangement  is  described. 
In  the  third  section,  the  model  development  is  presented. 
The  identification  scheme  is  discussed  in  the  fourth  section. 
Preliminary  studies  of  parameter  trends  are  presented  in  the 
fifth  section.  In  the  last  section,  remarks  are  collected  and 
presented. 

2.  Piezoelectric  microresonators 

The  microresonator  design  studied  in  this  work  consists  of  a 
clamped-clamped  microstructure  machined  out  of  a  composite 
wafer  [  I  ],  as  shown  in  figure  1 .  The  dynamic  characteristics  of 
the  resonators  provide  spectral  filtering  of  the  electrical  signals. 
In  order  to  produce  the  desired  piezoelectric  actuation  and 
sensing  capabilities,  the  device  is  fabricated  from  a  composite 
wafer  [14].  The  different  layers  of  one  style  utilizing  the 
piezoelectric  properties  of  lead  zirconate  titanate  (PZT)  are 
shown  in  figure  2.  The  suspended  structure  is  produced  by 
removing  a  portion  of  the  silicon  substrate  from  underneath 
the  oxide  layer.  This  silicon  dioxide  layer,  which  provides 
the  primary  stiffness  for  the  resonator,  offsets  the  neutral  axis 
of  the  structure  from  the  centerline  of  the  piezoelectric  film, 
and  the  two  platinum  layers  serve  as  electrodes  for  the  PZT. 
The  top  electrode  is  separated  to  create  input  and  output 
ports  for  the  microresonator.  Representative  dimensions  of 
geometry  corresponding  to  figure  2  are  presented  in  table  1 . 

In  a  second  type  of  microresonators  that  is  studied, 
aluminum  gallium  arsenide  (AlGaAs)  is  used  for  the 
piezoelectric  layer  [15].  For  this  style  of  resonator,  the  two 
electrode  layers  are  produced  by  doping  AlGaAs  with  silicon 
(AlGaAs:Si).  The  thickness  of  the  bottom  layer  of  AlGaAs:Si 
is  increased  in  this  design  to  provide  additional  stiffness  for  the 
resonator  and  eliminate  the  need  for  an  additional  layer,  such 
as  SiC>2  in  the  PZT  resonators.  Some  of  these  resonators  are 


Figure  3.  Structure  of  the  AlGaAs  microresonator. 
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Figure  4.  Experimental  arrangement  utilized  to  excite  and  monitor 
the  transverse  oscillations  of  a  microscale  resonator. 

Table  1.  Representative  PZT  resonator  dimensions. 


Dimension 

Symbol  Value 

Resonator  length 

/ 

200  /xm 

Width  of  all  layers 

b 

20  /xm 

Thickness,  Si(>2 

h\ 

1.06  /xm 

Thickness,  bottom  Pt 

b2 

135  nm 

Thickness,  PZT 

h 3 

530  nm 

Thickness,  top  Pt 

*4 

200  nm 

Table  2.  Representative  AlGaAs  resonator  dimensions. 


Dimension 

Symbol 

Value  (/xm) 

Resonator  length 

/ 

200 

Width  of  all  layers 

b 

10 

Thickness,  bottom  AlGaAs:Si 

A. 

2.0 

Thickness,  AlGaAs 

h2 

1.0 

Thickness,  top  AlGaAsiSi 

b 3 

0.5 

being  fabricated  to  be  thicker  than  the  PZT  resonators.  This 
results  in  a  stiffer  device  with  higher  resonance  frequencies, 
and  the  onset  of  nonlinear  behavior  in  this  device  occurs  at 
higher  excitation  levels.  This  design  also  helps  eliminate  some 
of  the  problems  that  can  arise  when  materials  with  different 
lattice  structures  are  placed  together  on  a  composite  wafer.  The 
arrangement  of  the  three  layers  of  the  AlGaAs  microresonators 
is  shown  in  figure  3.  One  set  of  dimension  values  for  this 
resonator  type  are  presented  in  table  2. 

The  arrangement  of  some  of  the  key  components  of 
the  experimental  setup  is  illustrated  in  figure  4.  With  the 
sample  fixed  in  an  RF  probe  station,  a  function  generator 
or  dynamic  signal  analyzer  is  used  to  apply  a  harmonic 
signal  while  the  response  is  monitored  by  using  laser 
interferometry.  After  initially  conducting  the  frequency  sweep 
in  a  quasi-static  fashion,  simulations  are  conducted  with 
the  identified  parameter  values  and  it  is  determined  that  a 
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Figure  5.  Frequency-response  data  for  a  200  /zm  PZT  resonator. 
Data  obtained  during  an  increasing  frequency  sweep  are  represented 
by  asterisks  and  data  obtained  during  a  decreasing  frequency  sweep 
are  represented  by  circles. 

swept-sine  signal  from  the  analyzer  can  be  used  to  obtain 
valid  results.  The  signal  produced  by  the  laser  vibrometer 
is  sent  back  to  the  analyzer  to  produce  a  frequency-response 
plot.  A  representative  frequency-response  data  set  is  shown  in 
figure  5,  where  the  response  amplitude  is  given  in  nanometers 
and  the  excitation  frequency  is  given  in  kilohertz.  Data 
represented  by  asterisks  4*’  correspond  to  an  increasing  (or 
forward)  frequency  sweep  and  data  represented  by  circles  ‘o’ 
correspond  to  a  decreasing  (or  backward)  frequency  sweep. 
Although  data  are  shown  for  both  forward  and  backward 
frequency  sweeps,  only  the  forward  sweep  data  are  used  in 
the  parametric  identification  scheme. 

3.  Resonator  model 

The  frequency-response  data  obtained  from  the  microscale 
resonators  clearly  reveal  nonlinear  characteristics.  Although 
the  dimensions  of  each  resonator  may  suggest  that  a  plate 
model  may  be  appropriate,  a  beam  representation  is  found  to 
adequately  model  the  dynamic  behavior  of  the  device.  To 
account  for  the  nonlinear  behavior  of  the  resonators  within 
the  model,  both  axial  stretching  and  nonlinear  curvature 
are  considered  [16].  By  applying  a  single-mode  Galerkin 
approximation  to  the  nonlinear  beam  model,  it  is  found  that 
the  contribution  of  the  nonlinear  curvature  to  the  coefficient 
of  the  resulting  cubic  term  is  more  than  three  orders  of 
magnitude  smaller  than  the  contribution  of  the  axial  stretching 
term.  For  this  reason,  the  nonlinear  curvature  is  considered 
to  be  negligible.  Because  of  this,  only  the  axial  stretching 
is  included  when  modeling  the  traverse  deflection  Wn  as  a 
function  of  axial  position  x  and  time  t. 

lPn(t)h(x)lxx  =  pAnWnM+cWnJ  +  (ElnWn'XX\xx  +  P0W„.*x 

j\wm.^<^Wn.xx.  (1) 

Within  (1),  the  parameters  and  variables  used  are  as 
follows:  mass  per  length  pAn ,  viscous  damping  coefficient 
c,  bending  stiffness  £/„,  compressive  axial  force  P0,  axial 
stiffness  EAm  and  axial  force  from  piezoelectric  layer  P„. 


Furthermore,  h(x)  is  a  discontinuous  function  that  describes 
the  separation  between  the  neutral  axis  and  the  applied  axial 
force.  This  separation  produces  the  distributed  moment  that 
excites  the  resonator  [7].  The  subscripts  n  and  m  are  used  to 
identify  the  sections  of  the  stepwise  axially  varying  geometry 
of  the  structure  and  M  is  the  total  number  of  sections.  The 
summation  is  required  to  calculate  the  total  axial  stretching 
within  all  of  the  sections.  Subscripts  following  a  comma  V 
indicate  partial  derivatives.  To  complete  the  model,  boundary 
conditions  and  compatibility  conditions  are  required.  The 
boundary  conditions  constrain  the  displacement  and  slope  of 
the  beam  at  each  end  to  be  zero.  Since  these  resonators  are 
made  up  of  three  different  segments,  the  displacement  profile 
is  separated  into  three  sections  that  require  continuity  at  the 
two  intersections.  In  order  to  accomplish  this  consistency, 
the  compatibility  conditions  require  that  the  displacement, 
slope,  moment  and  shear  be  equal  on  both  sides  of  each  of  the 
connections.  The  composite  structure  of  the  resonators  also 
requires  additional  calculations  to  obtain  averaged  properties 
such  as  the  bending  stiffness  and  the  axial  stiffness.  These 
values  are  calculated  based  upon  the  assumption  that  the 
effect  of  the  coupling  between  bending  and  stretching  is 
negligible.  With  the  nonlinear  partial  differential  equation, 
the  boundary  conditions,  the  compatibility  equations  and  the 
averaged  properties,  the  model  of  the  resonator  is  complete. 

In  order  to  facilitate  the  analysis  of  this  model,  the 
Galerkin  procedure  is  employed.  Noting  that  the  response 
range  of  interest  is  close  to  the  resonator’s  first  natural 
frequency,  a  single-mode  approximation  is  used.  Through 
this  procedure,  the  simplified  model  takes  the  form  of  a  forced 
Duffing  oscillator,  shown  as  (2)  where  z(t )  is  the  temporal 
amplitude.  The  coefficients  of  this  equation,  to  be  identified 
through  the  parametric  identification  scheme  are  the  following: 
the  modal  mass  m,  the  viscous  damping  coefficient  c,  the  linear 
stiffness  coefficient  k ,  the  nonlinear  stiffness  coefficient  a 3  and 
a  modal  force  parameter  F. 

With  this  model,  a  number  of  the  parameters  can 
be  predicted  by  using  the  resonator’s  geometry,  material 
properties  and  the  resonator’s  modeshape  function.  With 
either  a  target  frequency  or  an  independent  measurement  of 
the  axial  force  Pq ,  the  modeshape  can  be  determined  and 
used  to  calculate  the  values  of  m,  k  and  ot$.  An  approximate 
value  of  Pn ,  and  subsequently  F,  can  be  calculated  by  using 
the  piezoelectric  property  of  resonator  and  the  magnitude  of 
the  applied  sinusoidal  signal.  This  is  discussed  in  depth  in 
section  5.  With  the  aid  of  the  parametric  identification  scheme 
to  be  presented  in  section  4,  the  parameters  m,c,k,a 3,  F  and 
Pq  can  be  identified  from  experimental  frequency-response 
data  of  the  piezoelectric  microscale  resonators. 

mz{t)  +  cz(t )  +  kz(t)  +  of3  z3(/)  =  F  cos(ct >/).  (2) 

The  overdots  in  (2)  are  used  to  represent  time  derivatives. 
A  harmonic  excitation  is  assumed  in  arriving  at  (2).  The 
different  modal  coefficients  are  defined  by 


E,\[’  d*| 

(3) 

:Xl|/  J 

(4) 
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k  =  I  I  <t>n(x)[EI„<t>lv (X)  -  P0<t>:w]  dx  j 
n  =  l  • 

“3  =  {  f 

^(,)  =  E(/  " 0n(jc)  { ^[/>"<r)fi(jr)1} djr)  * 


(5) 


(6) 

(7) 


where  the  primes  in  (5)  and  (6)  indicate  spatial  derivatives  and 
<p„(x)  is  the  considered  modeshape  function.  Previous  studies 
have  shown  that  a  forced  Duffing  oscillator  with  a  hardening- 
type  nonlinearity  is  capable  of  producing  a  frequency-response 
curve  with  a  structure  similar  to  the  experimental  frequency- 
response  data  obtained  from  the  piezoelectric  microscale 
resonator  [17].  Assuming  the  nonlinearity,  damping  and 
forcing  to  be  weak  and  focusing  on  excitation  frequencies  close 
to  the  first  natural  frequency,  the  method  of  multiple  scales 
is  used  to  obtain  an  approximate  solution  to  the  nonlinear 
differential  equation  (2)  [18].  The  approximate  solution  for 
this  case  is  found  as  (8)  where  HOT  stands  for  higher  order 
terms  and  the  amplitude  and  phase  are  governed  by  (9)  and 
(10),  respectively. 


z(t)  =  a(t)  cos(a)t  -  y(t))  +  HOT  (8) 

fl(/)  =  -/za(0  +  *  sin(y(/))  (9) 

aO)y(t)  =  cra(t)  -  |cra3(/)  +  K  cos (y(t)).  (10) 

Periodic  responses  of  the  microresonator  correspond  to 
the  fixed  points  (a0.  Yo)  of  (9)  and  (10),  that  is,  a(t )  =  y(t)  = 
0.  The  fixed-point  equations  provide  the  frequency-response 
equation  (11),  which  shows  how  the  amplitude  of  the  periodic 
response  changes  with  respect  to  the  excitation  amplitude  and 
the  excitation  frequency.  The  parameter  values  at  which  the 
fixed  points  lose  stability  arc  given  by  the  critical  points’ 
equation  (12). 

[m2  +  [o  -  \aal)2]al  =  K1  (11) 

M2  +  (a  ~  l“flo)(°  “  laao)  =°-  (12) 

The  variables  within  (11)  and  (12)  are  related  to  the 
parameters  of  (2)  by  the  following  relations: 

con  =  Jk/m,  Q  =  o)/(On,  a  =  Q  —  1 , 

(13) 

[i  =  c/(2mcoo),  a  =  a3/k .  K  =  F/(2k). 

The  critical  points,  where  bifurcations  occur,  are  satisfied  by 
both  (11)  and  (12).  Analytical  curves  produced  by  using  these 
equations  are  displayed  in  figure  6. 

In  order  to  determine  the  spatial  function  needed  to 
calculate  the  coefficients  for  the  forced  Duffing  equation  from 
the  nonlinear  beam  model,  the  linearized  system  is  considered 
and  from  analysis  of  this  system  along  with  the  boundary  and 
compatibility  conditions,  the  first  natural  frequency  and  the 
associated  mode  shape  are  determined.  This  mode  shape  is 
used  as  the  spatial  function. 


Frequency 


Figure  6.  Thick  lines  are  used  to  show  the  analytically  predicted 
frequency-response  curve,  and  the  stable  and  unstable  response 
segments  are  shown  by  solid  lines  and  dashed  lines,  respectively. 
The  thin  lines  are  the  critical-point  curves,  which  are  independent  of 
the  excitation  level.  The  encircled  intersections  of  these  lines 
represent  the  critical  points  of  the  frequency-response  curves. 


Figure  7.  The  effects  of  increasing  parameter  values  on  the 
analytical  response  curve.  The  curve  depicted  by  using  solid  lines 
corresponds  to  the  nominal  case,  and  the  curve  illustrated  by  using 
dashed  lines  corresponds  to  the  effect  of  increasing  the  following: 
(a)  viscous  damping  coefficient,  ( b )  nonlinear  stiffness  coefficient, 
(c)  linear  stiffness  coefficient  and  ( d)  modal  force  parameter. 


4.  Identification  scheme 

With  the  aid  of  the  nonlinear  beam  model,  the  system 
parameters  are  examined  to  determine  how  they  affect 
the  structure  of  the  frequency-response  curve.  In  order  to 
determine  how  best  to  design  the  parametric  identification 
process,  the  damping  coefficient,  the  linear  stiffness 
coefficient,  the  nonlinear  stiffness  coefficient  and  the  modal 
force  parameter  are  examined.  The  influence  of  the  modal  mass 
is  not  examined  because  the  modal  mass  is  calculated  from 
the  nominal  geometry,  material  properties  and  approximated 
mode  shape  of  the  resonator.  Frequency-response  curves 
showing  the  effects  of  an  increase  in  each  of  these  parameters 
are  displayed  in  figures  1(a)  through  1(d).  As  expected 
and  shown  in  figure  1(a) ,  the  damping  coefficient  affects 
the  amplitude  of  the  peak  of  the  frequency-response  curve 
while  leaving  much  of  the  rest  of  the  curve  unchanged. 
The  nonlinear  stiffness  coefficient  influences  the  amount  by 
which  the  peak  of  the  frequency-response  curve  leans  away 
from  the  neutral  position  corresponding  to  the  first  natural 
frequency,  as  illustrated  in  figure  1(b).  Changing  the  value  of 
the  linear  stiffness  results  in  a  shifting  of  the  system’s  linear 
natural  frequency  and  this  causes  the  horizontal  position  of  the 
peak  of  the  frequency-response  curve  to  change,  as  shown  in 
figure  7(c).  The  modal  force  represents  the  general  magnitude 
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Figure  8.  A  flow  chart  depicting  an  overview  of  the  parametric 
identification  scheme  utilized  to  quantify  the  behavior  of  the 
microscale  resonators.  Starting  from  the  top  left  comer, 
experimental  data  and  initial  parameter  value  guesses  are  input  into 
the  identification  scheme.  The  three-part  identification  process  is 
employed  to  tune  the  parameter  values  by  comparing  the 
experimental  data  with  an  analytical  model.  These  three  steps  are 
repeated  until  the  desired  precision  is  attained.  The  resulting 
parameter  set  is  varied  and  reanalyzed  to  ensure  that  they 
correspond  to  the  global  optimum. 

of  the  force  being  applied  to  the  system  and  a  change  in  the 
value  of  this  parameter  affects  the  response  amplitude  of  the 
entire  frequency-response  curve  most  notably  at  the  peak;  this 
is  illustrated  in  figure  1(d). 

4.1.  Stages  of  identification 

Based  on  the  information  presented  in  figure  7,  the  parametric 
identification  scheme  is  developed  so  that  the  parameters 
with  the  most  influence  over  the  response  are  fit  before  the 
parameters  that  affect  the  response  the  least.  A  flow  chart 
depicting  an  overview  of  the  parametric  identification  scheme 
is  presented  in  figure  8.  First,  it  is  necessary  to  select  initial 
values  for  each  parameter  so  that  the  peak  of  the  analytical 
curve  is  located  within  the  same  frequency  range  as  the 
experimental  data  and  exhibits  a  similar  structure.  While 
the  selection  of  these  values  has  a  very  large  effect  on  the 
results  of  the  parametric  identification  process,  the  initial 
values  are  easily  determined  by  using  relations  obtained  from 
the  Galerkin  procedure  and  by  using  a  visual  comparison  of 
the  predicted  analytical  curve  and  the  experimental  data. 

Once  the  initial  values  are  selected,  a  three-stage 
parametric  identification  scheme  is  applied  to  determine  the 
optimal  combination  of  parameter  values  in  order  to  fit  the 
analytical  curve  to  the  experimental  data.  In  the  first  stage, 
the  frequency-response  equation  (11)  is  utilized  for  a  least- 
squares  curve  fitting  process.  The  parameters  determined  in 
this  section,  in  the  order  identified,  are  the  linear  stiffness 
coefficient,  the  nonlinear  stiffness  coefficient  and  the  modal 
force  parameter.  These  parameters  are  identified  by  tuning 
each  parameter  as  long  as  the  difference  between  the  analytical 
curve  and  the  entire  experimental  data  set  continues  to 
decrease.  The  parameters  are  tuned  both  positively  and 
negatively  for  a  number  of  different  adjustment  sizes. 

In  the  second  stage  of  the  identification  process, 
the  frequency-response  equation  (11)  as  well  as  the 
critical  points  equation  (12)  are  utilized  to  fit  the  upper 


bifurcation  point  by  tuning  the  viscous  damping  coefficient. 
As  previously  mentioned,  the  identification  procedure  is 
designed  specifically  to  study  nonlinear  dynamic  behavior 
that  experiences  jumps  in  the  frequency-response  data. 
Experimental  data  sets  that  do  not  experience  jumps  can  be 
studied  with  one  of  the  identification  schemes  previously 
discussed,  such  as  the  method  developed  by  Malatkar  and 
Nayfeh  ( 1 3].  Here,  a  similar  tuning  procedure  is  used  but  only 
the  analytically  predicted  amplitude  of  the  upper  bifurcation 
point  is  compared  with  the  amplitude  of  the  last  experimental 
data  point  before  the  jump.  Although  this  does  not  take  into 
account  any  difference  in  the  frequency,  additional  iterations  of 
the  identification  process  account  for  this  discrepancy  through 
tuning  of  the  linear  stiffness  coefficient.  While  both  the  modal 
force  parameter  and  the  viscous  damping  coefficient  affect  the 
amplitude  value  at  the  peak  of  the  frequency-response  curve, 
only  the  modal  force  parameter  significantly  affects  the  portion 
of  the  curve  away  from  the  peak.  As  a  result,  it  is  possible  to 
determine  a  unique  combination  of  these  two  parameters  for  a 
given  set  of  experimental  data. 

In  the  final  of  the  three  stages  within  the  identification 
scheme,  the  linearized  system  is  used  to  obtain  an  approximate 
mode  shape.  This  mode  shape  is  determined  by  tuning  the 
axial  force  within  the  linearized  system  to  synchronize  the 
first  resonance  frequency  of  the  model  with  the  frequency 
calculated  from  the  identified  parameter  values.  The  axial 
force  value  is  tuned  in  the  same  fashion  as  the  other  parameters. 
This  provides  an  approximation  of  the  microresonator’s  mode 
shape  for  this  natural  frequency  and  an  axial  force  value  that 
can  be  used  to  calculate  the  average  residual  stress  in  the 
resonator.  By  using  this  approximate  mode  shape,  the  modal 
mass  value  can  be  recalculated.  This  new  modal  mass  value  is 
then  used  in  the  next  iteration  of  the  identification  process  to  aid 
in  the  further  refinement  of  the  values  of  the  other  parameters. 

Upon  completion  of  these  three  stages  of  the  identification 
scheme,  the  parameter  values  are  examined  to  determine  if 
the  desired  level  of  precision  has  been  obtained.  If  the 
parameter  values  are  not  determined  to  sufficient  precision, 
the  identification  scheme  is  returned  to  the  first  stage  of  the 
three-stage  process  by  using  the  current  parameter  values  for 
the  initial  values.  Since  the  system  is  nonlinear  and  has  a 
multi  dimensional  parameter  space,  the  possibility  exists  that 
the  identification  scheme  may  converge  to  local  optima  and 
not  the  global  optimum.  To  avoid  selecting  local  optima,  after 
the  scheme  has  gone  through  a  sufficient  number  of  iterations 
to  produce  the  desired  level  of  precision,  each  of  the  key 
parameter  values,  c,k,ct}  and  F,  is  perturbed  both  positively 
and  negatively  to  produce  eight  different  sets  of  parameter 
values.  These  eight  sets  are  then  used  as  initial  values  and 
the  results  of  the  nine  cases  are  compared.  If  the  first  set  of 
parameter  values  is  found  to  have  the  best  values,  then  it  is 
considered  to  represent  the  ‘global  optimum’.  The  quality 
of  each  set  of  parameter  values  is  determined  by  calculating 
an  RMS  error  from  comparisons  of  the  resulting  analytical 
values  with  the  experimental  data.  In  the  event  that  the  best 
set  is  one  of  the  eight  perturbed  sets,  these  optimized  values 
are  perturbed  to  produce  eight  new  sets  and  the  same  process 
is  repeated. 


1597 


A  J  Dick  et  al 


Frequency  (kHz) 


Frequency  (kHz) 


Figure  9.  Nonlinear  frequency-response  curves:  (a)  a  200  /iw  PZT  microresonator  and  (b)  a  200  pm  AlGaAs  microresonator. 
Experimental  data  are  represented  by  asterisks  and  circles.  Analytical  curves  that  are  plotted  with  solid  and  dashed  lines  correspond  to 
parameters  obtained  with  the  parametric  identification  scheme  and  the  nonlinear  beam  model,  respectively. 


4.2.  Identified  parameters  for  PIT  and  AlGaAs 
microresonators 

The  same  experimental  arrangement  and  procedures  are 
employed  for  both  PZT  and  AlGaAs  microresonators.  In 
figure  9(a),  a  comparison  of  the  experimental  data  and 
frequency-response  curves  obtained  on  the  basis  of  the 
identified  parameter  values  and  beam  model  parameter  values 
is  shown  for  a  PZT  microresonator.  The  large  difference 
between  the  beam-model  response  and  the  identified  response 
is  reflected  by  the  parameter  values  and  this  is  discussed 
shortly.  For  the  PZT  microresonator,  a  minimal  RMS  error 
value  of  1.0943  nm  is  calculated  from  the  forward  sweep 
data  and  the  identified  parameter  values.  Although  only  the 
forward  sweep  data  are  used  to  determine  the  parameter  values, 
the  analytical  curve  plotted  with  identified  parameter  values 
in  figure  9(a)  shows  agreement  with  both  the  forward  sweep 
data  and  the  backward  sweep  data.  For  this  particular  data 
set,  the  RMS  error  calculated  for  the  backward  sweep  data 
has  an  even  smaller  value  of  1.0814  nm.  This  indicates 
that  although  the  backward  sweep  data  are  not  taken  into 
account,  the  identification  scheme  is  able  to  predict  the  correct 
frequency  response  in  the  lower  branch,  including  the  jump 
location. 

In  figures  9(6),  the  experimental  data  collected  from  a 
200  /xm  AlGaAs  resonator  are  presented  with  analytical  curves 
produced  with  the  parameter  values  identified  by  using  the 
authors’  identification  scheme  and  those  calculated  by  using 
the  nonlinear  beam  model.  The  difference  seen  between  the 
beam  model  response  and  the  identified  response  is  similar 
to  that  seen  for  the  PZT  device.  Once  again,  the  parametric 
identification  scheme  is  able  to  match  the  analytical  curve  to 
the  experimental  data  well.  For  the  data  shown  in  figure  9(6), 
an  extremely  small  RMS  error  value  of  0.5666  nm  is 
calculated. 

In  tables  3  and  4,  the  parameter  values  determined  by 
using  the  parametric  identification  scheme  and  numerical 
values  calculated  from  the  nonlinear  beam  model  are  shown 
for  the  PZT  microresonator  and  AlGaAs  microresonator, 
respectively.  The  second  column  of  each  table  contains 
parameters  calculated  by  using  the  identified  axial  force 
value.  The  parameters  in  the  third  column  correspond  to 
residual  stress  levels  determined  by  using  the  wafer  bow 


Table  3.  Comparison  of  parameter  values  for  the  PZT  resonator. 


Parameters 
and  units 

Parametric 

identification 

Beam  model 
Beam  model  (bow  stress) 

m  (kg  x 10-M) 

1.58 

1.58 

1.58 

c  (N  s  m-1  xlO-8) 

7.81 

- 

- 

k  (N  m_l) 

60.91 

60.88 

69.78 

«3  (Nm'3  xl0+l2) 

32.4 

3.54 

3.53 

P(Nx  lO”3) 

1.50 

1.50 

1.87 

(Do  (kHz) 

312.61 

312.55 

334.61 

GefTeclive 

397.15 

- 

- 

<W.ge  (MPa) 

43.56 

43.56 

54.29 

Table  4.  Comparison  of  parameter  values  for  the  AlGaAs  resonator. 


Parameters 
and  units 

Parametric 

identification 

model 

Beam 

Beam  model 
(bow  stress) 

m  (kg  xl0-M) 

1.13 

1.13 

1.12 

c(Ns  m"1  xlO-8) 

1.55 

- 

- 

J^Nm"1) 

54.02 

54.02 

45.83 

<*3  (Nm'3  xl0+l2) 

24.6 

2.84 

2.84 

P  (N  xlO-4) 

-2.68 

-2.68 

-5.88 

O)o  (kHz) 

347.58 

347.58 

321.42 

^effective 

1599.19 

- 

- 

Average  (MPa) 

-8.94 

-8.94 

-19.6 

measurement  method  (e.g.,  [19]).  From  comparisons  of  the 
identified  resonance  frequencies,  it  can  be  confirmed  that  the 
AlGaAs  microresonator  is  stiffer  than  the  PZT  microresonator. 
The  experimental  AlGaAs  data  also  show  lower  response 
amplitude  values  than  the  PZT  resonator  data  for  equivalent 
input  signals.  This  observation  is  consistent  with  the  fact 
that  the  AlGaAs  resonators  are  stiffer  and  experience  smaller 
transverse  deflections  and  that  the  piezoelectric  coupling  of 
the  AlGaAs  material  is  not  quite  as  strong  as  that  of  the 
PZT  material.  However,  as  seen  from  tables  3  and  4,  the 
AlGaAs  resonator  has  a  higher  effective  quality  factor  Q 
compared  to  the  PZT  resonator.  This  higher  (^-amplification 
for  the  AlGaAs  resonators  compensates  substantially  for  the 
increased  stiffness  and  reduced  piezoelectric  coupling. 

In  the  case  of  the  PZT  microresonator,  for  the  identified 
axial  force,  the  numerically  obtained  parameter  values  from 
the  beam  model  generally  agree  with  the  identified  values  with 
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the  largest  discrepancy  being  between  the  values  obtained 
for  the  nonlinear  stiffness  coefficient.  For  the  numerical 
values  calculated  with  the  residual  stress  from  the  wafer  bow 
measurements,  the  increased  stress/axial  force  values  affect 
the  linear  stiffness  coefficient,  and  subsequently,  the  first 
resonance  frequency. 

For  the  AlGaAs  microresonator,  there  is  also  a  difference 
between  the  residual  stress  levels  that  affect  the  linear 
stiffness  coefficient  and  frequency.  When  comparing  the 
parameter  values  within  table  4,  the  deviation  between  the 
identified  and  calculated  values  for  the  nonlinear  stiffness 
parameter  is  consistent  with  results  obtained  for  the  PZT 
resonators.  A  number  of  possible  sources  of  this  discrepancy 
are  investigated.  A  single-mode  Galerkin  approximation  is 
selected  because  the  frequencies  examined  are  in  close 
proximity  to  the  resonators  first  natural  frequency  and  the 
second  natural  frequency  has  a  value  of  more  than  twice  the 
value  of  the  first.  Studies  have  shown  that  the  use  of 
^he  domain  method  with  the  device’s  modeshapes  is  one  of 
the  more  successful  means  to  produce  a  reduced  order  model 
[20].  In  addition,  a  profile  assembled  from  multiple  single  point 
measurements  is  found  to  most  closely  match  the  first  linear 
modeshape.  When  conducting  the  perturbation  analysis,  a 
weak  nonlinearity  is  initially  assumed  and  found  to  match  well 
to  the  experimental  data.  Plotting  the  RMS  error  against  the 
excitation  amplitude  shows  that  low  excitation  levels  produce 
nearly  linear  behavior,  moderate  excitation  levels  produce 
weakly  nonlinear  behavior  and  higher  levels  of  excitation 
produce  stronger  nonlinear  behavior  that  diverges  from  the 
model  and  produces  higher  RMS  error  values.  The  curvature 
nonlinearity  is  not  included  in  the  nonlinear  beam  model  for 
the  reason  stated  in  section  3.  The  consistent  deviation  between 
the  nonlinear  beam  model  and  the  experimental  data  suggests 
that  there  may  be  an  additional  source  of  nonlinearity  that  must 
be  taken  into  account  in  the  nonlinear  beam  model  to  more 
accurately  represent  the  dynamic  behavior  of  the  considered 
microresonators. 

Overall,  as  demonstrated  here,  it  has  been  possible  to 
produce  a  successful  curve  fit  to  the  experimental  data  by 
using  the  parametric  identification  scheme  constructed  on  the 
basis  of  the  Duffing  oscillator. 

5.  Trends  of  identified  parameters 

Here,  results  of  preliminary  investigations  conducted  into  the 
variations  of  the  identified  parameters  with  respect  to  the 
selected  inputs  or  operating  conditions  are  reported.  By 
studying  the  parameter  trends  caused  by  changes  in  the  input 
signal  and  operating  conditions,  it  is  possible  to  determine 
how  to  improve  the  resonator  design  and  how  to  preprocess  the 
input  signal  to  improve  the  resonator  performance.  Different 
swept-sine  signals  with  increasing  amplitudes  and  the  addition 
of  various  dc  bias  levels  to  fixed  amplitude  sinusoidal  signals 
are  considered.  The  application  of  a  dc  voltage  has  been  used 
to  control  the  stiffness  of  a  MEMS  device  [21].  The  addition 
of  a  dc  bias  to  the  excitation  signal  is  also  used  to  account 
for  remanent  polarization  of  the  piezoelectric  material.  The 
parameters  examined  include  the  modal  force  parameter,  the 
linear  stiffness  parameter  and  the  nonlinear  stiffness  parameter. 


Figure  10.  Approximating  a  linear  force-voltage  relation.  The 
circles  ‘o’  represent  modal  force  values  for  different  data  sets  and 
the  solid  line  is  a  linear  curve  fit  to  these  data  points. 

5.1 .  Modal  force  and  piezoelectric  coefficient 

In  order  to  understand  how  changes  in  the  amplitude  of 
the  applied  signal  affect  the  modal  force  parameter,  the 
corresponding  term  in  the  nonlinear  beam  model  is  examined. 
The  excitation  results  from  the  distributed  moment  produced 
by  the  piezoelectric  material.  The  excitation  term  in  (1)  is 
separated  into  a  time-dependent  function  P„(t)  and  a  position- 
dependent  function  h(x).  Following  the  application  of  the 
Galerkin  procedure,  the  modal  force  function  F(t)  given  by 
(7)  is  produced.  After  carrying  out  integration  by  parts  twice 
with  respect  to  x,  (14)  is  obtained. 

F(0  =  P\{t)h4>\(xx).  (14) 

By  using  the  block  force  method  to  obtain  the  axial  force 
from  the  piezoelectric  material,  the  definition  of  the  modal 
excitation  term  can  be  further  expanded  to  include  material 
properties,  additional  device  geometry  and  the  applied  voltage. 
The  block  force  method  is  a  highly  approximate  method  that 
treats  a  piezoelectric  actuator  as  a  line  force  without  accounting 
for  any  spanwise  variation  of  stress  or  strain.  This  form  of  the 
modal  force  is  given  in 

F{t)  =  EMdu/DVinh^xi).  (15) 

This  equation  can  be  used  to  gain  addition  information 
about  the  device.  From  a  number  of  data  sets  collected  from 
the  same  device  for  increasing  drive  voltage  values,  a  complete 
set  of  parameter  values  can  be  identified.  A  plot  of  the  modal 
force  parameters  versus  the  drive  voltage  amplitude  yields  a 
near-linear  set  of  points  under  ‘ideal’  conditions.  An  example 
of  this  data  is  displayed  in  figure  10.  Rearranging  (15)  to 
solve  for  the  piezoelectric  coefficient,  du,  reveals  a  term  that 
directly  corresponds  to  the  slope  of  the  data  points.  Because 
this  term  consists  of  two  temporal  functions,  it  is  important  to 
note  that  the  slope  of  the  data  points  represents  the  ratio  of  the 
amplitudes  of  the  two  harmonic  functions  and  a  negative  sign 
must  be  included  to  account  for  the  phase  difference.  This 
form  of  the  equation  is  shown  in  (16),  where  the  actuator  area 
A  =  bh,  Fo  is  the  amplitude  of  F(t)  and  V0  is  the  amplitude 
of  V(t). 

<h\  =  -(F0/V0)(Ebh^(x]))-'.  (16) 

With  additional  information  about  the  resonator,  it  is 
possible  to  calculate  values  for  the  piezoelectric  coefficient. 
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Figure  11.  Trends  of  the  first,  linear  natural  frequency  for  different  PZT  microresonators.  Diamonds  correspond  to  100  /zm  resonators, 
squares  correspond  to  200  /zm  resonators  and  circles  correspond  to  400  /zm  resonators.  Identified  parameter  values  are  presented  in  (a)  and 
model  parameter  values  are  presented  in  ( b )  showing  how  the  model  is  capable  of  explaining  the  parameter  trends. 


Figure  12.  Trends  of  the  nonlinear  stiffness  coefficient  for  different  PZT  microresonators.  Identified  parameter  values  are  presented  in 
(a)  and  model  parameter  values  are  presented  in  ( b )  showing  how  the  model  is  partially  capable  of  explaining  the  parameter  trends. 


This  is  useful  since  material  properties  of  thin-film  materials 
often  differ  from  their  macroscale  counterparts  and  the 
properties  of  these  materials  can  be  dependent  on  the  various 
fabrication  procedures  to  which  they  are  exposed.  In  initial 
applications  of  the  identification  scheme,  the  determined 
piezoelectric  coefficient  values  range  from  —116  x  10-12^ 
to  -192  x  10-12^.  These  thin-film  material  values  are  of 
the  same  order  of  magnitude  as  the  piezoelectric  coefficient 
values  associated  with  bulk  ceramic  materials.  However, 
the  thin-film  material  piezoelectric  coefficients  are  smaller  in 
magnitude. 

5.2.  Axial  force,  linear  stiffness  and  nonlinear  stiffness 

The  addition  of  a  dc  bias  to  the  applied  signal  produces  a 
constant  axial  force  in  the  microresonator.  As  the  value  of 
the  dc  bias  increases,  the  peak  of  the  frequency-response 
curve  shifts  to  the  right.  Based  on  the  parametric  study 
previously  discussed,  the  axial  force  produced  by  the  addition 
of  the  dc  bias  affects  the  linear  stiffness  of  the  resonator. 
This  behavior  agrees  with  the  fundamental  understanding  of 
the  effects  of  an  axial  force  within  clamped-clamped  beam 
structures.  By  using  the  nonlinear  beam  model,  it  is  possible 
to  calculate  linear  natural  frequency  values  that  qualitatively 
and  quantitatively  agree  with  the  frequency  changes  observed 
in  the  experimental  data.  While  the  model  utilized  within 
this  work  does  not  include  the  effects  of  the  hysteresis  of 
the  piezoelectric  material,  a  basic  linear  approximation  is 
compared  with  identified  parameter  values  to  obtain  a  rough 
range  for  the  dc  bias  and  the  axial  force.  By  using  the  block 
force  method,  that  is,  axial  force  =  EA(^3i//)(dc  Bias), 


refined  values  of  this  voltage-force  ratio  can  be  used  to 
determine  the  piezoelectric  coefficient  d^\  for  the  devices. 
By  using  the  ratio  values  from  this  basic  model,  piezoelectric 
coefficient  values  ranging  from  —93.3  x  10“12^  to  -163  x 
10-12^  are  determined  for  the  PZT  microresonators. 
Although  this  is  merely  a  coarse  approximation,  these  values 
are  of  a  magnitude  common  for  the  piezoelectric  coefficients 
of  bulk  PZT  materials. 

To  qualitatively  compare  the  identified  parameters  with 
those  of  the  nonlinear  beam  model,  the  changes  in  the 
identified  first  natural  frequencies  are  shown  in  figure  1 1(a) 
for  a  range  of  dc  bias  levels  and  the  changes  in  the  first  natural 
frequencies  obtained  from  the  model  are  shown  in  figure  1 1  ( b ) 
for  a  range  of  axial  force  values.  Within  these  figures,  the 
diamonds  ‘O’  correspond  to  a  100  /zm  PZT  resonator,  the 
squares  *□’  represent  the  frequency  values  of  a  200  /zm  PZT 
resonator  and  the  changes  in  the  frequency  values  of  a  400  /zm 
PZT  resonator  are  represented  by  the  circles  ‘o’. 

Another  parameter  that  is  significantly  affected  by  the 
addition  of  a  dc  bias  is  the  nonlinear  stiffness  parameter. 
Increasing  the  level  of  dc  bias  causes  the  identified  nonlinear 
stiffness  parameter  values  to  decrease.  The  nonlinear  stiffness 
parameter  values  also  appear  to  be  influenced  by  the  hysteresis 
of  the  piezoelectric  material.  Again,  by  using  only  a  basic 
model,  a  range  of  axial  force  values  is  found  where  the  value 
of  the  nonlinear  stiffness  decreases  as  the  axial  force  increases. 
To  qualitatively  compare  the  identified  parameters  with  those 
of  the  nonlinear  beam  model,  the  identified  nonlinear  stiffness 
values  are  shown  in  figure  12(a)  for  a  range  of  dc  bias  levels 
and  the  nonlinear  stiffness  parameters  from  the  model  are 
shown  in  figure  12 (b)  for  a  range  of  axial  force  values,  both 
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normalized  with  respect  to  their  initial  parameter  values.  The 
axial  force  values  are  normalized  with  respect  to  the  buckling 
force  so  that  buckling  occurs  at  a  normalized  axial  force  value 
of  negative  one. 

In  this  section,  parameter  trends  with  respect  to  the 
resonator  inputs  have  been  examined  by  applying  the 
parametric  identification  scheme  to  a  series  of  data  sets.  These 
parameter  trends  included  modal  force  parameter  values  for 
various  drive  voltage  amplitudes  and  linear  and  nonlinear 
stiffness  variations  with  respect  to  dc  bias.  By  using  the 
nonlinear  beam  model,  it  has  been  possible  to  explain  the 
trends  of  each  of  the  parameters  and  produce  numerical  values 
that  qualitatively  agree  with  the  identified  values.  While 
further  work  will  be  necessary  to  develop  more  complete 
models  of  the  trends,  identifying  this  qualitative  agreement 
is  important  in  validating  the  selection  of  the  nonlinear  beam 
model  and  determining  the  course  of  future  work. 

6.  Concluding  remarks 

In  this  paper,  a  parametric  identification  scheme  is  developed 
with  the  capability  of  analyzing  nonlinear  systems  that  exhibit 
jumps  in  their  frequency-response  behavior.  In  addition  to 
its  ability  to  identify  the  parameter  values  corresponding  to  a 
nonlinear  model,  this  identification  scheme  can  also  provide 
additional  information  about  a  device  such  as  the  average 
residual  stress  level.  The  identification  scheme  is  applied 
to  frequency-response  data  collected  from  piezoelectric 
microscale  resonators  in  order  to  quantify  their  nonlinear 
behavior.  Methods  are  developed  to  determine  the  values 
of  an  equivalent  viscous  damping  coefficient,  a  linear  stiffness 
coefficient,  a  nonlinear  stiffness  coefficient,  a  modal  force 
parameter  as  well  as  the  level  of  axial  force  and  the  modal 
mass  of  the  resonator.  The  identification  scheme  has  been 
successfully  applied  to  data  obtained  from  both  PZT  and 
AlGaAs  microscale  resonators.  The  parameter  values  are  also 
compared  to  numerical  values  obtained  from  a  beam  model 
and  agreement  is  seen.  The  identification  scheme  is  applied  to 
multiple  data  sets  and  parameter  trends  have  also  been  studied. 

In  future  work,  it  is  conceivable  that  a  parametric 
identification  scheme  such  as  that  discussed  in  this  work  can 
be  used  to  tailor  the  characteristics  of  a  MEMS  array,  so  that 
phenomenon  such  as  localization  [22,  23]  can  be  engineered 
to  improve  the  performance  of  the  considered  devices. 
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